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ABSTRACT 


In an attempt to mitigate the expanding counter-UAV capabilities of adversary countries de¬ 
veloped in response to the United States’ increased reliance on these platforms, we apply a 
nondeterministic search pattern to a finite area searcher. By implementing a Levy distribu¬ 
tion on search leg lengths we analyze the trade-offs between efficiency and evasiveness of the 
searcher, comparing the expected time to target detection for a given set of Levy parameters 
to a probabilistic time to counter-targeting based on intelligence driven enemy capability. The 
culmination of this thesis is the development of a robust simulation tool, capable of model¬ 
ing various parameters on both searcher and search area, the output of which is a quantifiable 
estimate on the probability of mission success. 


V 



THIS PAGE INTENTIONALLY LEET BLANK 


VI 



Table of Contents 


1 Introduction 1 

1.1 Background. 1 

1.2 Research Objective. 2 

1.3 Literature Review. 2 

1.4 Scope, Limitations, and Assumptions. 9 

1.5 Contributions of the Thesis. 10 

1.6 Organization of the Thesis. 10 

2 Mathematical Formulation 13 

2.1 Background. 13 

2.2 The Search Area and Searching Unit. 17 

2.3 The Nondeterministic Search Pattern. 21 

3 Modeling and Simulation 33 

3.1 Simulation Model Assumptions. 33 

3.2 Developing the Levy Distribution. 34 

3.3 Restriction of the Searcher to a Finite Area. 36 

3.4 Algorithm for Searcher Movement. 39 

3.5 Graphical Display. 43 

4 Results, Analysis, and Experimentation 51 

4.1 Probability of Mission Success. 51 

4.2 Determining Coverage Ratio. 52 

4.3 Incorporating Counter-Targeting. 55 

4.4 Providing Decision Support. 59 

4.5 Field Experimentation. 64 

vii 























4.6 Bayesian Update and Looping Function. 67 

5 Conclusions and Future Avenues 73 

5.1 Conclusion. 73 

5.2 Avenues for Future Work. 74 

5.3 Operational Impact and Recommendations. 75 

A Summary of Variables and Symbols 77 

Appendices 77 

B Software and Hardware 81 

C Waypoint Data Processing for Field Experiments 83 

C.l Waypoint Generation. 83 

C.2 Display of Waypoint Data. 83 

D Calculations for Discrete Sample Spacing 87 

List of References 90 

Initial Distribution List 99 


viii 









List of Figures 


Figure 1.1 Levy probability density plot . 5 

Figure 2.1 Motivational search scenario. 14 

Figure 2.2 (a) An example of perfect non-overlapping complete search, commonly 

referred to as the “lawnmower” or “spiral-out” pattern, (b) The resultant 
linear increase to total area coverage exhibited by the “spiral-out” search 
pattern. From [41]. 14 

Figure 2.3 (a) An example of the continuous random search, commonly referred to 

as the “confetti search.” (b) The resultant exponential increase to asymp¬ 
totically complete coverage of the total area exhibited by the “confetti 
search” pattern. From [42]. 15 

Figure 2.4 Efficiency vs. Evasiveness. 16 

Eigure 2.5 A boundary zone is established within the circular search area. 18 

Eigure 2.6 Sensor field-of-view determines footprint. 20 

Eigure 2.7 Simulation parameters and variables . 24 

Eigure 2.8 Dubins curve configurations. 26 

Eigure 2.9 Heading configurations for Dubins calculations, derived from [34] ... 26 

Eigure 2.10 Decision table for determining Dubins path, from [34]. 27 

Eigure 2.11 Illustration showing the variables required in the calculation of the Du¬ 

bins distance. This example shows the geometry for a RSL configuration, 
but the process is similar for LSR . 29 

Eigure 2.12 Illustration showing the configuration of pivot points at the initial and 
final waypoints of a search leg. The location of the pivot point is at 
a distance rtum from the waypoint, offset n radians in the appropriate 
direction. 30 
















Figure 3.1 Functional block diagram of simulation model. 33 

Figure 3.2 Matlab generated Levy probability density plots. 36 

Figure 3.3 Histogram of leg lengths used in look-ahead Levy search. 38 

Figure 3.4 Standard Levy CDF vs. look-ahead Levy CDF. 39 

Figure 3.5 Effects of discrete sampling on coverage. 39 

Figure 3.6 Graphical display of simulation model. 43 

Figure 4.1 A single run of the simulation model shows that during the early stages 

of the search, a coverage rate similar to that of the continuous ideal 
searcher is achieved. However, as more area is uncovered, the random¬ 
ness of the pattern causes the searcher to double back over areas previ¬ 
ously searched, resulting in an logarithmic shape to the curve. 53 

Figure 4.2 Plot of the mean coverage results from the n = 1000 simulation runs 
using the parameters in Table ??. The fitted regression curve is plotted 
for comparison. 55 

Figure 4.3 (a) Plot showing the effect on mission success of an increase in Rsearch, 

given scaling parameter, c, and Tiqck- (b) Plot showing the effect on 
mission success of an increase in scaling parameter, c, given Rsearch and 
Tlock- Changes in Rsearch have a greater effect on the probability of 
mission success. 57 

Figure 4.4 (a) Plot of probability of mission success versus search area radius for a 

given Tlock- (b) Plot of the probability of mission success versus Tiqck 
for a given search radius. 59 

Figure 4.5 A pairs plot of the data is generated in R to assist in quantifying the 
relationship between the search area size, target lock time, and the prob¬ 
ability of mission success. 60 

Figure 4.6 Histogram on time to target detection, Tq, resulting from n = 10,000 
simulation runs without enemy counter-targeting capability. The dark 
line function represents an exponential regression fit on the distribution. 61 

Figure 4.7 A variation of the “lawnmower” search discussed in Section ??, the 
“spiral-out” search begins at the origin and conducts straight-leg, pre¬ 
dictable, non-overlapping sweeps toward the search boundary edge. . . 62 


X 














Figure 4.8 This summary plot of the cumulative distribution functions for varying 
values of target lock time (n = 1000) can be used to determine a measure 
of the probability of mission success. 64 

Figure 4.9 Procerus Unicorn UAV [48] used in multi-UAV look-ahead Levy search 

demonstration conducted at 13-2 JIFX in February, 2013. 65 

Figure 4.10 Google Earth overlay of major Levy waypoints generated by LookA- 
headLevyFlight.m algorithm. The starting airfield, search center, and 
search boundary are annotated. 66 

Figure 4.11 Google Earth overlay of the actual flight path of a single UAV searcher 

conducting the look-ahead Levy search. 67 


Eigure 4.12 A magnified section of the flight path traveled by the UAV compared to 
the waypoints generated through the use of Dubins curves in LookA- 
headLevy.m. The simulation model created for this thesis is demonstra¬ 


bly an adequate representation of the flight physics of a small UAV. . . 68 

Eigure 4.13 Google Earth overlay of the actual flight path of a single UAV searcher 

conducting the look-ahead Levy search. 69 


Eigure 4.14 (a) Heat map showing the distribution of probability of target presence 

employing a perfect sensor Bayesian search over the area, (b) The prob¬ 
ability density of larger search area sections is used to bias the movement 
of the searcher upon expiration of a loop timer.. 70 


Eigure D. 1 This graph shows the effects of spacing between the discretized looks on 
the runtime of the simulation. As the spacing increases, the simulation 
runtime decreases exponentially. Choosing a discretized spacing equal 
to the search radius provides 96% of the coverage attained through a 
continuous search for a fraction of the runtime cost. 87 

Eigure D.2 This illustration shows the geometry involved in calculating the sensor 
overlap observed during a discretized search with spacing less than twice 
the sensor radius. 87 


XI 











THIS PAGE INTENTIONALLY LEET BLANK 



List of Tables 


Table 2.1 Variable associates with sensor operations. 18 

Table 2.2 Mid-size UAV operating parameters. 19 

Table 2.3 Variables associated with sensor operations. 22 

Table 2.4 Variables associated with searcher characteristics. 23 

Table 2.5 Variables associated with the calculation of Dubins curves and their path 

lengths. 28 

Table 4.1 Input parameters to simulation runs used to evaluate coverage rate. ... 53 

Table 4.2 Regression output from Excel Data Analysis tool. The low jp-value and 

high R-squared values indicate the goodness of fit to the coverage rate. . 54 

Table 4.3 Input parameters to simulation runs measuring probability of mission suc¬ 
cess against an anti-UAV asset with a Tiqck capability of 90 seconds. . 56 

Table 4.4 Results of n = 1000 simulation runs against an anti-UAV asset with a 
Tlock of 90 seconds. Mission success is achieved upon target detection, 
with mission failure occurring if the searcher is counter-targeted. 57 


Table 4.5 Sample of the results of n = 100 simulation runs against an anti-UAV as¬ 
set with variations on Tiqck and Rsearch for constant Levy distribution pa¬ 
rameters. Mission success is achieved upon target detection, with mission 
failure occurring if the searcher is counter-targeted. Parameters leading 


to a higher probability of mission success than mission failure are bold¬ 
faced. 58 

Table 4.6 Times to target detection utilizing a deterministic “spiral-out” search, de¬ 
pendent on distance of the target from the beginning of the search. ... 63 

Table 4.7 Input parameters for simulation runs to assess Td and Tc . 63 












Table 4.8 Searcher parameters established for real-world flight during multi-UAV 

JIFX demonstration. 66 

Table 4.9 Inputs to the simulation incorporating a Bayesian update and looping 

function. 69 

Table 4.10 Comparison of performance between the standard Look-ahead Levy search 

and a look-ahead Levy search implementing the looping function and 
Bayesian update. The looping function reduces the average time to tar¬ 
get detection, but the forced override of the length distribution sometimes 
leads to exceeding the critical length, resulting in counter-targeting of the 
searcher, and overall mission failure. 70 

Table A.l Simulation variables set by the user . 77 

Table A.2 Variables which are populated throughout the simulation with information 

useful to the analysis of search performance. 80 


XIV 








List of Acronyms and Abbreviations 


A2AD Anti-access/Area Denial 
AOU Area of Uncertainty 

ARSENL Advanced Robotic Systems Engineering Laboratory 
CDF Cumulative Distribution Function 

CRUSER Consortium of Robotics and Unmanned Systems Education and Research 

DHS Department of Homeland Security 

FOC Farthest On Circle 

JIFX Joint Inter-agency Field Exploration 

MANPADS Man-portable Air Defense System 

NPS Naval Postgraduate School 

PDF Probability Distribution Function 

SLS Sea Level Standard 

UAV Unmanned Aerial Vehicle 

UAS Unmanned Aerial System 

USG United States Government 


XV 




THIS PAGE INTENTIONALLY LEET BLANK 


XVI 



Q 


Executive Summary 




Unmanned Aerial Vehicles (UAVs) have become a mainstay of modem iday military operations, 
such as surveillance and reconnaissance missions. Adversaries of the United States are fully 
aware of this shift and are developing weapons and training to counter these unmanned assets. 
Many of the counter-UAV weapons under development and testing require that a target lock be 
maintained on the UAV-iQr some minimum amount of time. By randomizing the flight pattern, 
e.g., when conductingsKRViissions, and limiting the time the UAV travels on any single flight 
leg, we can minimize the vulnerability of these assets. To accomplish this randomization we 
employ a Levy distribution function to determine the length of each search leg, while changes 
in searcher heading are drawn from a uniform distribution. We model realistic flight limitations 
using Dubins curves, which define the minimum distance path between two points of differ¬ 
ent heading orientation given the minimum turn radius capability of the searcher. Regression 
analysis of simulated search times is used to derive the expected coverage rate. We define a met¬ 
ric, the probability of mission success, comprised of a time to target detection by the searcher 
and a time to counter-targeting of the searcher by the adversary. A Bayesian update scheme 
is applied to the search to incorporate imperfections in sensor performance, along with a loop¬ 
ing search function. Should no target detection occur within a specified amount of time, the 
searcher’s travel will be biased toward the area of highest target probability density. The culmi¬ 
nation of this thesis is the development of a simulation model for analysis on the employment of 
nondeterministic search patterns as a means to mitigate counter-targeting and counterdetection 
threats. 
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CHAPTER 1; 
Introduction 


1.1 Background 

Unmanned aerial vehicles (UAVs) have become the mainstay of modern-day intelligence gath¬ 
ering for military and government agencies worldwide. These assets provide valuable real-time 
situational awareness while minimizing several of the risks associated with personnel in poten¬ 
tially adversarial territories. A variety of mathematical models for searching such areas have 
historically been applied to provide decision support for these information-gathering missions. 
However, many search models and associated measures of performance address only the effec¬ 
tiveness of the search and deal solely with the constraints of the searcher’s characteristics, such 
as limited endurance or limited sensor performance. 

Instead, the mission performance may require that the employed search model address addi¬ 
tional concerns that are relevant to realistic mission scenarios, such as the need to operate in 
unfriendly or clandestine environments and avoid counter-detection and/or counter-targeting. 
The ability to encompass mission-relevant metrics in a holistic mathematical model, rather than 
to construct individual models for different components of a complex mission, offers significant 
benefits to decision makers who must assess the relevant trade-offs. 

For example, although maximal area coverage is essential, as nations develop more advanced 
air tracking and defense systems, the need for stealth by either platform or employment of 
randomized flight patterns that incorporate unpredictability (to avoid detection and targeting 
by the adversary) also increases, requiring a balance between search effectiveness and searcher 
evasiveness. Threats from emerging technologies such as directed energy weapons [1,2], as 
well as conventional ones such as shoulder-launched rockets, pose potential danger to aerial 
assets conducting surveillance or other stand-off missions. 

In this study, we analyze the employment of a nondeterministic search pattern over a finite area 
using randomly distributed leg lengths and heading angles. We are interested in the probability 
of mission success, i.e., the ability of the searcher to discover a target located in the finite area 
prior to being counter-targeted by the adversary. 
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1.2 Research Objective 

The focus of this thesis is the development of an analytic model and associated computer simu¬ 
lation with which to analyze the area coverage and target detection capabilities of an unmanned 
searcher conducting randomized search over a finite area. The effectiveness of such a search is 
measured against the known efficiency of a continuous sweeping search by deriving equations 
representing the probability of detecting a target uniformly distributed within the finite area, as 
well as the expected time elapsed prior to the detection. Utilizing this estimated time to target 
detection we apply assumptions regarding adversary targeting capabilities to formulate a prob¬ 
ability of mission success, e.g., the probability that the time required for the searcher to detect 
the target is less than the time required for an enemy to counter-target the searcher. 


1.3 Literature Review 

The field of robotic coverage and search has an extensive history, including substantial efforts 
to optimize the coverage patterns for different objective functions such as minimum time or 
coverage path, minimal overlap, increased robustness to localization errors, etc. [3-5]. The role 
of randomized coverage and search patterns has also been explored, with studies ranging from 
optimal random paths to traverse all nodes of a finite network, to random walks over infinitely 
two-dimensional area [6-10]. Specifically, several studies have been conducted over the years 
in an effort to determine the proper distribution of search leg lengths resulting in an efficient 
random search, where efficiency is measured by the searcher’s ability to cover a large area in a 
finite period of time [11,12]. The demonstration of this type of efficient random search in nature 
has been the focus of multi-disciplinary efforts with animal behavior modeling, e.g., [13-17]. 

Viswanathan et al. recognized this efficiency in the natural foraging habits of animals, deter¬ 
mining that several species exhibit Levy characteristics in their search patterns [18,19]. Given 
that the livelihood of these species depend on efficient search techniques for survival, the mim¬ 
icking of these habits should translate to an efficient (if not optimal) search pattern for most 
object search scenarios. Viswanathan et al. analyzed the search leg distances recorded dur¬ 
ing a three-month observation of albatross (five albatross conducting nineteen foraging jour¬ 
neys) [20]. Unlike what might be observed through demonstration of a random Gaussian or 
Raleigh distribution, no well-defined variance was observed among the search leg distances. 
Plots of the recorded data indicated a long-tailed power-law distribution of lengths correspond¬ 
ing to Levy flight motions [18]. In particular, the observed distribution of flight lengths was 
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characterized by 


( 1 . 1 ) 


p{ti)~h " 

where ti is the time between landings, and ji represents the fractal dimension of the long-tailed 
power-distribution with 1 < /i < 3. The distribution converges to a Gaussian for /i > 3 and does 
not correspond to a normalizable density function for /i < 1. Histogram plots of the nineteen 
individual foraging trips suggested an average /i ~ 2. 

Further works by Viswanathan et al. assert that there are several animal foraging patterns ob¬ 
served in nature (albatross, honeybees, aphids, drosophila, etc.) which mimic a Levy flight 
with /i ~ 2. The authors credit these multiple observances of the Levy distribution in nature to 
the fractal properties of the Levy search that lead to a large dispersion from the origin over a 
finite time period. This characteristic is particularly beneficial for swarms or groups of N Levy 
walkers versus N Brownian walkers [21]. The Levy walkers are able to cover a more widely 
spread area and diffuse more quickly, thus limiting competition between individual searchers. 
Conversely, Brownian searchers tend to remain clustered together for longer periods of time. In 
summary, it is shown that due to the low probability of revisiting a previously searched area, 
the Levy flight possesses an advantage over Brownian motion when the target density is sparse 
and the targets are randomly distributed, but makes little difference in small areas of high target 
density [18]. 

The authors’ extended analysis of the Levy search differentiates between destructive and non¬ 
destructive search methods, the former being one in which an area is precluded from future 
searches once it has been visited, the latter allowing for multiple visits to any given site. A 
set of simulations is conducted in which the searcher is attributed a fixed search radius repre¬ 
sentative of an animal’s visual range, conducting a continuous search while moving along the 
Levy path [19]. If a target is detected within this visual range, the searcher can change direc¬ 
tion and move directly to the target, eliminating the target from the search area. If no target is 
detected within the visual range of the searcher, or if successful target elimination is achieved, 
a random heading from a uniform distribution and a random search length from the power-tail 
Levy distribution are drawn and the searcher moves along this new path until a target is ac¬ 
quired or the end of that search leg is reached. Through comparison of the simulation results 
to the foraging data collected on honeybees, deer, and wandering albatross, the authors affirm 
the observance of /i ~ 2 for any search dimension where there exists no a priori knowledge of 
target location [18,19]. 
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Raposo et al. expand upon the work of Viswanathan et al. by adding a parameter, t, representing 
a regeneration time during which a visited site is not available for revisit. The authors show that 
for any power-law tailed distribution (Equation 1.1), with the optimal exponent restricted to 
1 < Hopt < 2, by relating the velocity of the forager, v, and the recovery time of the site, T, such 
that z = X’ where A is the average distance (based on a spatial Poisson distribution of target 
location) between targets, the optimal search strategy is defined by Hopt = 1+ (jx)- 

As T approaches 0 (non-destructive limit), the best search strategy is observed with jiopt ~ 2. 
This coincides with the conclusions of Viswanathan et al. As T approaches infinity (destructive 
limit), the best search strategy is observed using fiopt —)■ 1, suggesting upper and lower bounds 
on the mean flight length for a Levy search pattern [22]. 

Calitoiu and Milici demonstrate that the Levy distribution exhibits characteristics advantageous 
to groups or swarms of animals foraging over large areas [23]. The authors describe the ben¬ 
efits of the Levy flight over the Brownian walk as being derived from the fact that the search 
leg distances associated with Brownian motion are scaled to a set magnitude defined by the 
searcher, where as the Levy flight lengths are a function of target spacing. This is observed in 
nature through a species’ ability to adjust to changing environmental conditions and availabil¬ 
ity of food/prey. The authors employ the Levy flight pattern to simulate destructive and non¬ 
destructive searches using non-cooperative agents; searchers working independently to locate 
targets in a common area with the constraint that two agents cannot occupy the same vicinity at 
the same time. The authors draw a specific distinction between the utilization of a Levy flight 
vice a Levy walk for the simulation. The term “Levy flight” is used to describe a pattern in 
which the search is conducted only at the ends of each leg; searchers use a “move, stop and 
search” technique. In the Levy walk, the searchers are continually able to detect a target using 
a “move and search, stop” technique. Lor both the destructive and non-destructive simulations, 
a shift parameter of /i = 2 is used. The Levy distribution searches are compared to searches 
utilizing a step size equal to unity, consistent with Brownian motion. In both destructive and 
non-destructive simulations the Levy searchers traverse distances farther from the origin. How¬ 
ever, because of the non-continuous search characteristics exhibited during the Levy flight, the 
Brownian searchers are able to search more total area. The results of these simulations infer that 
a Levy distribution of search leg lengths results in a dispersion from the origin that is beneficial 
in areas of low target density. The employment of a Levy walk vice Levy flight would further 
increase the effectiveness of this search. 

Each of the aforementioned works show that a Levy distribution of search leg lengths often 
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results in a more rapid dispersion from the origin when eompared to Gaussian-based Brownian 
motion, leading to a more expanded area of eoverage and less loeal eompetition between entities 
searehing as a group. It is this dynamie aspeet of the Levy distribution that makes it an intriguing 
eandidate for use in randomization of a seareh pattern over a finite area. 

While widely studied with regard to natural behaviors of animal speeies, the Levy distribution 
is also the foeus of multiple roboties-based seareh strategies, being applied to seareh leg lengths 
of autonomous searehers attempting to loeated targets in a minimal duration of time [24,25]. 
Lenagh and Dasgupta eonduet simulations utilizing what they deseribe as linear Levy and 
looped Levy seareh patterns (deseribed shortly) employing multiple searehers to deteet a single 
target loeated within an unbound seareh area [26]. The linear Levy seareh employs seareh leg 
distanees drawn from the elosed-form Levy distribution of Equation 1.2 and illustrated in Fig¬ 
ure 1.1. The seareher travels these Levy lengths eondueting a eontinuous seareh for the target, 
eonsistent with the Levy walk deseribed previously. 

rr 



Levy length, I. 


Figure 1.1: Probability distribution for the closed-form Levy with shift parameter /r = 0. The smaller 
the scaling parameter, c, the sharper the probability peak and greater chance of returning shorter lengths. 


In this elosed-form Levy distribution L(/y;/i,c), the variable /i represents the shift or loeation 
parameter for the distribution of lengths for leg Ij, defining the minimum length returned. The c 
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represents the scaling parameter controlling the height of the curve, and thus the probability of 
any particular length being returned. Figure 1.1 shows that for smaller values of c the probabil¬ 
ity peak of the distribution is higher, such that the probability of obtaining a value close to the 
parameter ji is greater [26]. The authors define the change in orientation of the searcher only to 
occur at angles greater than | from the current heading, thus the initial heading of the searcher is 
random, drawn from the uniform distribution, U (0,2;r), with subsequent turn directions defined 
as f + 0, where 6 ~ t/(0, ;r). The looped Levy search utilizes the same distributions; however, 
if contact is not gained/held on any target after a finite time period, the searcher will return to 
its origin (if no target information is retained) or to the area of last known target location (if 
contact was held but subsequently lost) and recommence the search in a new random direction. 
The searcher is equipped with a perfect sensor, such that if the target is within a distance less 
than the sensor range from the searcher, the target is located. These simulations are conducted 
in three phases, the first being search for a stationary target that has relocated to a spot some 
distance from its last known location. The remaining phases involve localizing (phase two) and 
tracking (phase three) of mobile targets. Each of these simulations is conducted for varying 
values of parameters /i and c. The authors conclude that little difference exists between the 
performance of the looped and linear Levy models while searching for stationary targets. The 
functionality of the looped Levy to return to the last known target position limits the time the 
searcher travels a path of zero target density, resulting in increased performance over the linear 
Levy when searching for moving targets. Throughout the tracking phase the looped Levy search 
demonstrates a timelier re-acquisition of a moving target after contact is lost, and for all simu¬ 
lations, faster target detection occurs when the value of the shift parameter /i is set to equal the 
initial target distance from the searcher (though this requires a priori knowledge of the target’s 
location) [26]. This distribution function representing the Levy walk is the foundation of the 
non-deterministic search developed in this thesis. The target locations are randomly selected 
and unknown by the searcher, thus dictating the need for further analysis regarding the optimal 
shift (/i) and scaling (c) parameter values. 

Mantegna and Stanley explore the significance of Levy movements within bounded physical 
systems, such as turbulent fluids or single molecules embedded in a solid, in which there exists 
an unavoidable cutoff in distance traveled [27]. This research provides insight as the the effects 
of conducting Levy-based searches within a bounded area of operations. The authors refer to 
the random movement within systems of this type as truncated Levy flights and derive a mean 
and variance for this truncated Levy distribution for which the mover is physically limited in 
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distances it can move. Mantegna and Stanley conclude that, unlike the standard Levy, the dis¬ 
tribution of this truncated Levy does converge to a Gaussian distribution (i.e., exhibits behavior 
prescribed by the Central Limit Theorem) when the sample size exceeds 10^ to 10^. 

As referenced, there have been several studies addressing the optimization of random search 
patterns, i.e., a series of random straight line distances linked by random changes in direction. 
The majority of these studies establish evidence that the rapid dispersion characteristic of the 
Levy distribution provides the advantageous opportunity to achieve maximum area coverage 
with minimal revisit time. These previous works provide the initial parameters for our model 
and simulation. This thesis builds upon these studies by utilizing the closed-form Levy distri¬ 
bution (Equation 1.2) to simulate and analyze the coverage ratio, that is, the area searched as a 
function of time of flight, and probability of detection for a Levy walk pattern over a bounded 
area. To model the realism of UAV movement and sensor capabilities, several other works are 
reviewed. 

In 1957, Lester Dubins proved that the shortest distance between two points in two-dimensional 
space, each defined by a coordinate and heading, can be expressed as a single path configura¬ 
tion of the set D = {CSC,CCC}, where C represents a curve of minimum turn radius and S 
represents a straight line segment [28]. These can be specified further such that the short¬ 
est path between the two points must consist of some path configuration D, where D G D = 
{RSLj RSRj LSLj LSRj RLRj LRL}, with R and L representing right and left turns, respectively. 
By calculating the length of each of the three path segments, the optimal path, i.e., that with 
the minimum overall distance, can be determined. This path is known today as “Dubins path.” 
This work is extended in the development and analysis of Reeds car, a Dubins vehicle allowing 
for both forward and backward motion, and the Dubins airplane, adding altitude as a functional 
parameter in the determination of the Dubins path [29-32]. Parlangeli et al. further the work on 
the Dubins model by introducing a path comprising various waypoints, each of which must be 
visited, in order, by a mover bound by limited curvature [33]. 

Shkel and Lumelsky simplify the mathematics behind Dubins curves by proposing a method 
which allows for the selection of the shortest Dubins path from the set D without having to 
conduct the cumbersome segment length calculations performed by Dubins [34]. Using four 
quadrants akin to that of a polar coordinate system, the authors distinguish classes of initial and 
final heading configurations, realizing that not every path in D is feasible for every configuration 
of headings. Further, these classes are divided into equivalency groups, which account for 
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equal but opposite flight paths of orthogonal pairs of heading configurations. For example, 
a searcher at an initial point with heading ^ moving to a final point with heading 0 would 
conduct a RSL flight path. Had the initial heading been — ^ the searcher would travel the same 
distance but opposite path, or LSR. The authors define sixteen classes arranged in a four-by- 
four matrix of quadrants, with each row representing the quadrant of the initial heading and 
each column representing the quadrant of the final heading. The intersection of the rows and 
columns contain information regarding the appropriate Dubins path, eliminating several of the 
previously required calculations. All calculations conducted by Shkel and Lumesky assume a 
unit turn radius in order to simplify the results; thus the scalability of their equations must be 
analyzed in order to be applied to this thesis, in which a non-unit or variable turn radius can be 
set by the user of the simulation. 

The operational characteristics of the sensor utilized by the searcher necessarily affects the per¬ 
formance of the search strategy regardless of pattern employed. It is essential for operational 
considerations to account for possible imperfections in sensor performance in order to obtain 
a full appreciation for any search pattern efficiency. Several works address the significance of 
these sensor imperfections, i.e., false and missed detections [35,36]. Kress et al. established a 
discretized search grid with each cell assigned a prior probability of containing a target and no 
cell capable of housing more than one target [37]. They simulate a search utilizing probability 
thresholds in order to determine if a cell is clear or contains a target. A Bayesian update process, 
incorporating the probability of false alarms and missed detections, is used to determine poste¬ 
rior probabilities of target location. The simulated UAV searchers used in this previous work are 
non-adaptive, meaning that they can not adjust the predetermined flight path to accommodate 
for the updated probabilities. To account for this, the searchers are allotted a specific number of 
looks, each of which is used only if probability of a target being located in the visited cell does 
not exceed the threshold for positive detection or clearance of the cell. Chung and Burdick also 
incorporate a Bayesian update scheme into a discretized search pattern, applying a “decision¬ 
making” process to searcher movement. The discrete Bayesian update process is combined 
with a look-ahead search algorithm to bias the searcher movement toward a region within the 
discretized area presenting the highest belief of target presence [38]. The incorporation of this 
“decision making” process and Bayesian update of the search grid are utilized in this thesis to 
implement a looped Levy search, directing the searcher to an area of highest target probability 
if no target is found within a finite time period. 
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1.4 Scope, Limitations, and Assumptions 


The purpose of this thesis is the mathematieal derivation and subsequent eomputer simulation 
of a nondeterministie seareh algorithm that maximizes the area eovered during the seareher’s 
flight while limiting the probability of eounter targeting. The seope of this thesis is limited to 
answering the following questions. Can we: 

• determine an optimal distribution of random flight leg distanees and turn angles whieh 
improves upon the eoverage ratio aehieved during a eontinuous random seareh, while at 
the same time redueing the single leg eounter-targeting exposure of the seareher? 

• develop a tangible and effeetive means of measuring counter-targeting avoidance? 

• develop a scalable nondeterministie search algorithm capable of implementation via any 
unmanned system (whether aerial [UAV], ground [UGV], surface [USV] or subsurface 
[UUV])? 

• derive a single specific adjustment parameter that can scale the model for any of the above 
search platforms? 

• define and achieve an acceptable coverage ratio and probability of detection observable 
over multiple simulated scenarios, and then use these observed characteristics of the 
model to derive a probability of achieving a given coverage ratio or probability of de¬ 
tection for various scenarios? 

For the purposes of this research, it is assumed that the searchers are “on station” at the onset 
of the search pattern (i.e., take-off and landing of the UAVs are not accounted for). The model 
addresses the realistic flight limitations of a UAV by limiting the searcher to travel along Dubins 
paths, with all heading changes limited to a minimum turn radius. Environmental effects such 
as wind/current are assumed to have no effect on the flight trajectory of the searcher. Other than 
directional changes, the motion characteristics of the searcher will remain fixed (i.e., no speed 
or altitude changes) throughout the entirety of the search pattern. It is further assumed that the 
search sensors are mounted on some form of gimbal device, such that changes in vehicle list 
and trim that would otherwise shift the sensor footprint need not be modeled. With respect to 
the search area, it is assumed that no obstacles or terrain obstructions exist that would interfere 
with the path nor sensor of the searcher. When simulating multiple searchers, each unit is able 
to maintain constant and perfect communications with all other units and the central controlling 
station. Finally, it is assumed that at most one target exists per each discretized cell of the search 
area. 
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1.5 Contributions of the Thesis 

The research accomplished during the development of this thesis is beneficial to the United 
States military and other government agencies operating unmanned vehicles in support of in¬ 
telligence, surveillance, and reconnaissance missions. The main contribution of this paper is 
the development of a probabilistic model by which to analyze the trade offs present when using 
an autonomous searcher to find a stationary target in minimal time while avoiding counter¬ 
targeting by the adversary. We define a novel performance measure which captures this balance 
between search effectiveness and evasiveness of the searcher, which can be used as a basis for 
assessing search mission performance in bounded adversarial environments. We utilize a Levy 
distribution in the generation of nondeterministic search paths, and highlight the insights arising 
from the analysis of the search and counter-targeting performances. 

Further, this thesis research represents the author’s selection as a recipient of the Space and 
Naval Warfare Systems Center, Pacific (SSC Pacific) NPS Student Fellowship Program, which 
enables interactions with SSC Pacific laboratory researchers, and also provides substantial 
equipment and travel funds in support of enhancing the student’s research. 

Additionally, a portion of the thesis work presented herein, describing the use of the Levy- 
distributed search leg lengths with an instantaneously turning searcher, has been accepted via 
peer-review for publication and presentation at the 2013 IEEE International Conference on 
Robotics and Automation (ICRA) in Karlsruhe, Germany [39]. This conference represents a 
flagship venue for dissemination of advanced research in the international robotics and automa¬ 
tion communities. 


1.6 Organization of the Thesis 

The remainder of this thesis is organized as follows. In Chapter 2, we include a detailed defini¬ 
tion of all parameters and assumptions utilized in the model and simulation. The derivation of 
the analytic model to determine coverage ratio with respect to time over a finite search area is 
addressed. Chapter 3 details the development and implementation of the computer simulation 
based on the analytic model, to include the incorporation of Dubins paths, Bayesian probabil¬ 
ity updates, and the looping functionality of the search. Chapter 4 analyzes the results of the 
simulation in comparison to expected values calculated in Chapter 2. Sensitivity analysis is 
performed on each of the adjustable input parameters of the model, in particular the shift and 
scaling parameters of the Eevy distribution and the size of the search area and sensor footprint. 
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In Chapter 5, the final conclusions are revealed, to include recommendations on the optimized 
values of the model parameters required to achieve the desired balance of searcher coverage 
and stealth. In this chapter we also address recommendations for follow-on studies with regard 
to the Levy search pattern. 
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CHAPTER 2: 
Mathematical Formulation 


This chapter describes the mathematical mechanisms and theory behind the Levy-based search 
pattern that is the cornerstone of this research. The reasons for the establishment of a finite 
search space to include the operational significance of the shape and size of the area are ad¬ 
dressed. A detailed description of the searcher movement and sensor characteristics based on 
parameters of existing real-world unmanned aerial craft is included. This chapter also describes 
the implementation of a Bayesian updating scheme in which prior and posterior probabilities 
of target presence are calculated based on sensor operating characteristics, contributing to the 
efficiency of target detection. 

2.1 Background 

This research can best be summarized through description of a search scenario encountered 
often during military operations. Assume an intelligence report is received by an operational 
unit indicating that a high value target is, or recently was observed, conducting operations at 
the location indicated in Figure 2.1. It is desired that the unit conduct a UAV search to confirm 
presence of the target for facilitation of a possible strike mission. The unit is not to violate 
the airspace of Country X, as indicated by the shaded area in the upper left. To further com¬ 
plicate the situation. Country X has established several anti-UAV assets along its borders and 
in the vicinity of the last known target position. Successful target detection requires that the 
unmanned searcher execute an efficient search while maintaining enough randomness of flight 
to avoid being counter-targeted and destroyed. A common measure of efficiency in search is 
the coverage rate of the searcher, or how much of the area can be searched in a given amount of 
time. In general, this rate, F{t), is a function of the speed and sensor footprint of the searcher, 

where v is the searcher velocity and w is the search sweep width of the sensor [40]. The most 
efficient search utilizes a sweeping pattern of area coverage as shown in Figure 2.2. Assuming 
the searcher is equipped with a perfect sensor, meaning that any target located within the max¬ 
imum range, rsensor, of the sensor is detected, the sweep patterns are such that no gap (missed 
coverage) and no overlap (excess coverage) exists between sweeps. For this searcher with a 
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- Restricted Air Space 

Aliied UAV 
Anti-UAV weapon 


Figure 2.1: The purpose of this research is to develop a nondeterministic search pattern capable of 
efficiently locating a target while minimizing vulnerability to enemy attack. 


perfect sensor conducting a continuous non-overlapping search of the total area, the coverage 
ratio with respect to time, F{t), is defined as 

F{t) = (2.2) 

^search 

where t is time diWAAsearch is the total area to be searched. Thus, at any time t, F{t) represents 
the percentage of Asearch searched. This continuous ideal search is an important benchmark 
in that no segment of the search area is missed and no segment is searched more than once, 
resulting in the efficient linear rise to total coverage shown in Figure 2.2. 


Coverage Faaot for Coninuous Complete Coverage 




Figure 2.2: (a) An example of perfect non-overlapping complete search, commonly referred to as the 
“lawnmower” or “spiral-ouf’ pattern, (b) The resultant linear increase to total area coverage exhibited 
by the “spiral-out” search pattern. From [41]. 


14 


































In contrast to the continuous ideal search is the continuous random search pattern. Random 
movement over the search area coupled with discrete “looks” results in an inefficient search 
pattern albeit asymptotically complete. This search is often likened to dropping equally sized 
pieces of confetti onto a surface, one at a time, until the entire surface is covered [40]. As with 
the ideal search, the time to complete coverage is a time-based function of the searcher speed 
and sweep width. However, because the movement over the search area is random, the potential 
exists for revisit of any segment in the area, resulting in a logarithmic shape to the coverage 
ratio over time. As more of the area is searched, the probability of a random movement to an 
area not previously covered decreases and the coverage exhibits a slow and asymptotic rise with 
“diminishing returns” to complete coverage. The equation for this continuous random coverage 
ratio is defined as: 

_ vwt 

F(t) = 1 —e ^search ^ (2.3) 

with vwt determining the size of the “confetti piece” uncovered with each “look” [40]. Illustra¬ 
tions of this search pattern and the resulting exponential shape of the coverage ratio due to the 
potential for revisit are shown in Figure 2.3. 



Figure 2.3: (a) An example of the continuous random search, commonly referred to as the “confetti 
search.” (b) The resultant exponential increase to asymptotically complete coverage of the total area 
exhibited by the “confetti search” pattern. From [42]. 


The continuous ideal search pattern demonstrates great exploration potential, efficiently cov¬ 
ering the search area in a minimal amount of time. However, the employment of a perfect 
sensor required to achieve coverage which exhibits no overlap and no excess is unrealistic when 
translating to real-world applications. The continuous random search shows great evasiveness 
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capability due to the unpredictable nature of the random movement, but the requirement to in¬ 
stantaneously move about the search area is again not transferable to real-world applications. 
The intent of this research is to design a search algorithm that promotes a balance between the 
efficient exploration characteristics of the ideal search and the evasion potential of the searcher 
as exhibited by the random pattern. This hybrid pattern can then be executed to perform the 
mission objective: maximize the probability of the UAV searcher detecting the target prior to 
being counter-targeted and destroyed by the anti-UAV assets (Figure 2.4). 



Figure 2.4: The goal of this research is to define a reafistic and reafizabfe search pattern promoting the 
efficacy of the ideaf search (dotted) white minimizing the evasiveness of the searcher through incorpora¬ 
tion of random movement. 


We develop a model by which to quantify this expected probability of mission success for the 
unmanned unit conducting a nondeterministic area search within some bounded search region. 
We apply tangible measures of performance (MOP) to this search mission in order to quantify 
the probability of successfully detecting a target. We define a time to target detection {Td) as the 
time elapsed from the start of the search until positive detection of the target. This parameter is a 
function of the searcher’s movement and sensor capabilities, the size of the search area, and the 
random distributions of heading angle and transit leg length used to govern the trajectory. An 
estimated time to counter-targeting {Tc) is derived from models of intelligence assessments of 
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adversarial technologies operating in and around the search area. This intelligence assessment 
consists of the estimated time required for an anti-UAV weapons system to achieve a target 
lock and subsequent firing solution on the searcher. Tc is the time elapsed from the start of the 
search until the time the searcher is destroyed due to traveling a predictable path, such as straight 
flight legs, for a time exceeding the assessed enemy capability. Given that Tc is a function of 
the counter-targeting capacity of the adversary, and operationally beyond our control, we treat 
it as a random variable, establishing distributions on counter-target times specific to various 
levels of capability. The total probability of mission success can then be derived from these two 
parameters, where mission performance is quantifiable, defined as the expected time to target 
detection being less than the time for counter-targeting of the searcher: 

P{To<Tc) Measure of Performance (2.4) 


2.2 The Search Area and Searching Unit 

2.2.1 Establishment of a Finite Search Area 

This research utilizes a finite search area designed to represent a circular area of uncertainty 
(AOU). An AOU is used in contact management situations when uncertain or time-late intelli¬ 
gence regarding a target’s position is obtained, or when contact on a known target is lost. The 
shape of the AOU is often circular or elliptical, with the boundaries defined by factors ranging 
from confidence in the intelligence source to the attributes of the target itself. Often the radius 
of this AOU is a function of the time late condition of the intelligence report and the maximum 
transit speed of the target, with the boundary referred to as the farthest-on circle (FOC). In other 
words, the size of the AOU is dependent on how far the target could have possibly traveled in 
the time since the location was reported. In many cases, it is probable that the target is not 
transiting at maximum capable speed on a single straight-line heading, and thus, the established 
AOU is a reasonable and realistic upper bound on the search area. The size and shape of the 
AOU may also be a product of the confidence in and fidelity of the intelligence report. Other 
factors such as how long, i.e., time late, it takes to position a searcher in the area must also be 
taken into account and can often increase the size of the AOU. 

In many scenarios, the air space or ground territory where the search is conducted poses sen¬ 
sitive territorial issues, requiring the establishment of a strict search boundary. To model this 
operationally relevant constraint, we establish a boundary zone, shown in Figure 2.5. The width 
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of this boundary zone is equivalent to twice the turn radius (2rturn) of the searcher. The maxi¬ 
mum straight distance between the initial and final waypoints, (xo^yo) and (xf,yf), respectively, 
of any search leg, will not exceed this artificial boundary, ensuring that the searcher will never 
violate the actual area boundary. This constraint still enables the searcher to move over the 
entirety of the search area. The parameters associated with the search area are summarized in 
Table 2.1. 



Figure 2.5: A boundary radius is established within the circular search area, ensuring that the searcher 
does not violate the actual area boundary. 


Variable 

Definition 

Units 

^search 

^search 

^boundary 

total search area 

radius of the circular search area 

effective radius of the search area, outside of which the UAV 
must be turning 

[meters^] 

[meters] 

[meters] 

Table 2.1: Variables associated with sensor operations 


2.2.2 The Searcher 

Searcher Movement 

The searcher modeled throughout this research is intended to represent UAV technology cur¬ 
rently in service, in particular, mid-sized platforms utilized by Special Forces and other forward 
deployed personnel. These systems are often rugged, low-cost UAVs capable of short range 
search and small payload delivery missions. Table 2.2 summarizes common operating parame¬ 
ters exhibited by these mid-size platforms [43]. 
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Parameter 

Capability 

Stall Speed 

18 m/s (40 MPH) 

Maximum Speed 

36 m/s (80 MPH) 

Cruise Speed 

22 m/s (50 MPH) 

Operational Maximum Range (still air) 

320 km (200 miles) 

Operational Endurance 

4 hours 


Table 2.2: General operating parameters of a mid-sized UAV, derived from [43,44] 

Sensor Footprint 

It is assumed for the purposes of this research that the UAV is equipped with a sensor package 
capable of distinguishing between the target and the environment in which it is operating. We 
assume a circular sensor footprint for ease of calculations. The field of view (FOV), that is, 
the extent of the target area visible during any instant in time, is assumed to be 30 degrees. 
This is derived from the works and reports of the U.S. Army Research Center and the Defense 
Technical Information Center and considered to be the average capability of small and mid¬ 
sized UAVs [44]. The “on the ground” footprint of the sensor is determined using a height of eye 
formula (Equation 2.5) based on the searcher altitude {alt) to calculate the sensor radius, rsensor, 
and thus, area of the circular footprint (Equation 2.6), as shown in Eigure 2.6 and Table 2.3. This 
sensor range, rsensor, when extended along both sides perpendicular to the searcher’s trajectory, 
determines the sweep width, w, referred to in the search theory examples of section 2.1. 

f sensor — ^ ifll^ ) tan(E OV ) 

^sensor = sensor 


(2.5) 

( 2 . 6 ) 


To simplify the promulgation of searcher and target position, as well as calculations on distance 
and probability of target presence, the search area is discretized into one meter by one meter 
square cells. Assuming no intelligence suggesting otherwise, each of these discretized cells has 
an equal opportunity of containing the target, assumed to be no larger than the one meter by one 
meter cell. As the UAV searcher moves over the search area, the sensor footprint uncovers these 
cells. Eor the purposes of this research, it is assumed that if the circular footprint of the searcher 
overlaps any part of a grid cell, the cell is considered searched. This assumption is plausible 
due to the large scale of the search area (generally of the magnitude 10^ square meters) as 
compared to the area of the searcher footprint (on the order of 10^ square meters). The possible 
over-accounting for the search of grid cells not completely covered by the sensor is therefore 
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FOV 



Figure 2.6: The search area uncovered during any discrete instant in time, Asensor, is a function of the 
sensor radius, Vsensor, as determined hy the sensor field of view, FOV, and the altitude of the searcher, 
alt. 


negligible, yet greatly simplifies the ealeulations involved when eomputationally determining 
whieh eells have been searehed. 


Modeling of Sensor Imperfections 

The assumption of a perfeet sensor is suffieient when determining the effeetiveness of the 
seareher movement pattern, but is not a feasible assumption for modeling real-world seareh 
dynamies. In assuming a perfeet “eookie eutter” sensor, any target that lies within the maxi¬ 
mum range of the sensor is seen. Due to adjustability in sensor gains and possible limitations 
in proeessing power for these mobile sensor systems, there exists, in reality, the tendeney to 
generate false eontaet reports or missed target deteetions. We denote the probability that a sen¬ 
sor indieates the deteetion of a target when one is not present, generally referred to as a false 
alarm, by the designation a. The probability that the sensor fails to indieate the presenee of a 
target when one is present, eommonly referred to as a missed deteetion, is denoted /3. A perfeet 
sensor is eharaeterized by both a and /3 equal to zero. 

Throughout this researeh, it is assumed that the target exists within the eonfines of the eireular 
seareh area; thus, the sum of the probabilities that the target exists in eaeh grid eell must equal 
one, i.e., 

1^ = 1, (2.7) 

i=\ 

where pi is the probability that the target is present in eell i and n is the total number of grid 
eells in the seareh area. In the absenee of any intelligenee information indieating that the target 
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is in a certain part of the search area, the probability that the target is in each cell is initialized 
as 

Pi=-: Vz = (2.8) 

n 

to represent a uniform probability. To derive the probability that the target exists in each individ¬ 
ual grid cell, these probabilities of false alarms and missed detections can be applied according 
to a Bayesian updating scheme. Bayesian updates use the a priori probability that the target 
exists in a cell, pu and apply the a and sensor characteristics to determine an a posteriori 
probability the target exists. If a searcher detects a target in a given cell, the following update is 
applied to that cell’s target probability [38]: 


P 


ftarget 
I present 


‘target 

detected” 


(l-i8)pz 

{a){\-pi) + {\-^)Pi' 


(2.9) 


As more detections occur in that cell, pi increases per Equation 2.9, until it exceeds a predefined 
threshold probability, at which point it can be deduced that the target is actually located 

in cell z. In a similar manner. Equation 2.10 utilizes the lack of detection as the searcher passes 
over a given cell to update pi, with the probability of target existence decreasing for that cell 
until falling below 


P 


j target 
I present 


“target not 
detected” 


mPi 

{I - a){l - Pi) + (15) Pi' 


( 2 . 10 ) 


Eor the purposes of this research, the target is assumed stationary; therefore, the overall prob¬ 
ability that the target exists in the search area must remain at 100 percent. Thus, after the 
Bayesian update is accomplished for each individual cell, that change in probability must be 
distributed to all other cells in search area. In this way, the probability of target existence in 
cells where detections occur increases and in cells with no detections decreases, regardless of 
whether those cells have been physically searched yet. Table 2.3 summarizes the variables 
associated with the operations of the sensor. 


2.3 The Nondeterministic Search Pattern 

2.3.1 Development of the Levy Search 

The Eevy-based search utilized throughout this thesis is a trajectory of the searcher comprising 
piecewise linear segments of Eevy-distributed lengths, as exhibited in numerous natural sys¬ 
tems [18,20,21]. Such a model offers an initial mathematical basis for exploring the balance 
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Variable 

Definition 

Units 

FOV 

field of view of the UAV onboard sensor 

[degrees] 

alt 

altitude of the searcher 

[meters] 

r sensor 

sensor operating range of detection 

[meters] 

^sensor 

circular area of sensor footprint 

[meters^] 

a 

sensor probability of false alarm 

[none] 

j8 

sensor probability of missed detection 

[none] 

Pi 

current probability that a target is in a specific cell i 

[none] 

con firm 

Pi 

probability above which a grid cell is confirmed to contain 
a target 

[none] 

deny 

P i 

probability below which a grid cell is assumed clear of a 
target 

[none] 


Table 2.3: Variables associated with sensor operations, as shown in Figure 2.6 


between ideal search (for improved search effectiveness) and random search (for improved un¬ 
predictability in motion). We refer to the relationship of these attributes as the efficacy versus 
the evasiveness of the searcher, respectively. Upon initialization of the Levy search, the searcher 
is located at the starting position, (xo^yo), the Cartesian equivalent of the latitude and longitude 
representing the center of an area of uncertainty (AOU) in which it is desired to conduct the 
search. A random heading of the j* iteration, fj, is drawn from a uniform distribution on turn 
angles from 0 to In. The individual length of the y* flight leg, Ij, is drawn from the Levy 
distribution shown in Figure 1.1 and formulated (repeating Equation 1.2) as 

rv 

Recall that the scaling parameter, c, affects the uniformity of the distribution on lengths whereas 
the shift parameter, /i, limits the shortest length sampled from the distribution. The heading 
angle of each subsequent search leg, Yj+i^ is determined by adding a random heading change 
6, again drawn from the uniform distribution from 0 to In, to the current searcher heading fj, 
as in Equation 2.12. 

Yj+i = Yj + 0, with 0 ~ t/(0,2;r). (2.12) 
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The searcher travels from {xo,yo) to {xf,yf), the point located length Ij from the initial point 
and calculated as 


Xf = Xo + ljCos{Yj), (2.13) 

yf = yo + ljSin{Yj). (2.14) 

A summary of these parameters is shown in Table 2.4 and illustrated by Figure 2.7. 


Variable 

Definition 

Units 

Ij 

length of single flight leg j, drawn from 

[meters] 

Yj 

initial heading of searcher on leg j, drawn from U (0,2;r) 

[radians] 

e 

change of heading from legy to leg j+i 

[radians] 

Yj+i 

heading on leg y -l- 1, where 7y+i = Yj + 0 

[radians] 

V 

speed of searcher 

[meters/sec] 

r turn 

turn radius of the searcher 

[meters] 

^max 

{xo,yo) 

maximum endurance of the searcher 

Cartesian position of searcher at beginning of search leg y 
Cartesian position of searcher at end of search leg y 

[seconds] 


Table 2.4: Variables associated with searcher characteristics 


2.3.2 Bounded Levy Walk 

The Levy walk has been shown to perform well while conducting a search over an infinitely 
large area due to the rapid dispersion from the search origin, resulting in larger area coverage 
over time [26]. In real-world operations, the search platform involved will have limited en¬ 
durance capability and will generally be constrained to a finite area of operations defined by the 
AOU or FOC referenced in Section 2.2.1. To maximize the potential for the searcher to locate 
the target within these constraints, a reasonable estimate of target location must be known and 
bounds placed on this operational search area. This research explores two methods to confine 
the search pattern to the vicinity of the established AOU. The main intent of these methods 
is preserve searcher endurance by limiting the amount of excess coverage (i.e., coverage ap¬ 
plied to areas outside of the AOU), but these bounding behaviors can also ensure the searcher 
is contained within the flight space cleared for the search, e.g., preventing a breach of territorial 
boundaries in contested areas. We introduce these methods here, and address them in detail in 
the following chapter. 
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Figure 2.7: Illustration of the relevant parameters in determining the next search leg of the simulated 
Levy flight. 


Turn-around 

We refer to the first of the boundary enforeement methods as the “tum-around” method. As 
depleted in the basie description of the search pattern, a random leg length Ij is drawn from the 
Levy distribution function and applied to the flight path on a random heading a. After moving 
that leg of the search path, the “tum-around” algorithm compares the current searcher location 
to the center of the AOU. If it is determined that the searcher is beyond the bounds of the AOU, 
a bias is placed on the random headings, 6, ensuring the searcher moves back toward the AOU 
center. Once it is determined the searcher is back within the boundary of the search area the 
bias is lifted and the random search continues. This method allows the searcher to temporarily 
exit the search area, resulting in excess coverage and waste of precious endurance capacity and 
mission performance, but preserves the underlying Levy distribution for leg lengths. 


Look-ahead 

The second version is referred to as the “look-ahead” method. Through implementation of this 
method, a leg length and heading are drawn from the associated distributions and a calculation 
is performed, prior to travel, to determine the endpoint of the new leg. If the new heading and 
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angle combination places the searcher outside the bounds of the area, that Ij is discarded and 
both a new length and heading angle are drawn, repeatedly so until the sampling results in a 
terminal position that is in-bounds. 


2.3.3 Dubins Levy Walk 

The Levy-based search algorithm uses an instantaneous change of direction to transition from 
leg to leg. This is unrealistic when attempting to model the physical characteristics of a fixed 
wing UAV, which must constantly be moving forward at speed greater than the designed stall 
speed. We must account for the minimum turn radius capable by the UAV during a heading 
change in order to obtain a realistic portrayal of the actual flight path and area covered during the 
search. We model this realistic flight trajectory through the implementation of Dubins curves. 

In 1957, Lester Dubins developed a system of equations to determine the optimal shortest path 
of a constant mover between any two waypoints, regardless of the heading configuration at each 
waypoint [28]. He derived that the path between these points could be broken down into three 
component parts: a starting turn, a straight-line path or second turn, and a finishing turn. The 
configuration of the three parts is based on the distance between the two points, the minimal 
turn radius of the mover, and the initial and final heading of the mover at the waypoints. By 
proving that each path can be broken down into these three path components, a total of six 
different paths exist. The set of these paths is known as Dubins set, which we annotate as D. A 
visualization of these configurations is shown in Figure 2.8. Each of the turns is denoted by an 
R for right or L for left and the straight-line path is denoted as an S, such that the set of behaviors 
is: 

D = {RSR,LSL,RSL,LSR,RLR,LRL}. (2.15) 


In 2001, Shkel and Lumelsky analyzed the complex trigonometry conducted by Dubins and 
were able to categorize all Dubins paths between two points into Table 2.10, comprising six¬ 
teen configurations of initial and final mover headings [34]. By imposing a constraint on the 
mover that the distance between any two waypoints must be greater than four times the mini¬ 
mum turn radius of the mover, the authors ensure the turn circles never overlap (Figure 2.8[e] 
and [f]), eliminating the LRL and RLR sequences from the set of possible Dubins curve config¬ 
urations. Depending on the quadrant associated with the initial and final angle headings (see 
Figure 2.9),there only exists up to four possible Dubins paths between the points, and in several 
cases, only one or two paths. These simplifications greatly reduce the number of calculations 
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(a) (b) 


0<J GH3 

(c) (d) 




Figure 2.8: Graphical representation of the possible comhinations of component parts of Duhins curves, 
illustrating trajectory configurations: (a) RSR (h) LSL (c) RSL (d) LSR (e) RLR (f) LRL. 


required to find the optimal (shortest) path between the two waypoints. 
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Figure 2.9: Shkel and Lumelsky determined that hy limiting the distance between the two points to 
no less than four times the turn radius and determining the quadrant of the initial {\xo,yo\, ct) and final 
{[xf,yf],l5) heading angles, the set of Duhins curves possible for each move is greatly simplified, derived 
from [34]. 
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,Tina|^Quadrar^ 
Initial Quadrant 

1 

2 

3 

4 

1 

RSL 

if S|2^ 0 then RSR 
if S|2^ 0 Ihcn RSL 

if S|j< 0 then RSR 
if S|j> 0 ihcn LSR 

if S[^> 0 then LSR 
if S|4 > 0 then RSL 

else RSR 

2 

if $21 < 0 ihcn LSL 
if $21 > 0 then RSL 

if SJ2 ** 0 thCT LSL 
if Sjj > 0 then RSL 
if S22 0 then RSR 

if S^, > 0 then RSL 

RSR 

if S24 0 then RSR 
if S24 > 0 then RSL 

3 

if S31 < 0 then LSL 
if S31 > 0 then LSR 

LSL 

if SJ3 < 0 then RSR 
if SJ3 > 0 then LSR 
if S33 < 0 then LSL 
if Si, > 0 then LSR 

ifSj^ < 0 then RSR 
ifSj4 > 0 then LSR 

4 

if S41 > 0 then RSL 
if sji > 0 then LSR 
else LSL 

if S42 ^ h Ihcn LSL 
if S42 ^ 0 Ihcn RSL 

if S43< 0 then LSL 
if S43> 0 then LSR 

LSR 


Figure 2.10: Decision table developed by Shkel and Lumelsky for determining the shortest path between 
two points of known heading. From [34]. 


The initial and final heading eonfigurations are determined using Figure 2.9. The distanee 
ealeulations in [34] are eonducted for eaeh of the possible Dubins eurves refereneed from Table 
2.10. The shortest of these distanees represents the length, J,of a nominal Dubins path of 
eonfiguration D G D, and is the sum of eaeh of the three eomponent parts; s^, the are length 
along the first turn, 5®, the length along the straight-line path, and s^, the are length along the 
seeond turn [34], sueh that: 

d = + + DgD. (2.16) 

In this researeh, we further simplify the ealeulations of optimal path eonfiguration by ensuring 
the heading of the seareher at the initial waypoint is always in the direetion of the seeond 
waypoint, sueh that the relative heading between the two points is zero. This eonstraint reduees 
the determination of optimal Dubins path to be dependent only on the final heading of the 
seareher, eliminating the need for a RSR or LSL turn eonfiguration. In summary, the eonstraints 
on heading and distanee between the two points allows for the elimination of the RSR, LSL, 
RLR, and LRL eonfigurations from the set of possible nominal paths. The seareher utilized in 
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this research only needs to conduct a LSR or RSL Dubins curve to achieve the nominal path 
between the two waypoints. 

Formulation of Dubins Curves 

This section provides a brief review of the geometric derivation of Dubins curves and calcu¬ 
lation of their associated path lengths, so as to facilitate their implementation in the proposed 
simulation framework. A summary of all variables used in the calculation of Dubins curve 
lengths is shown in Table 2.5.As mentioned previously, only the RSL and LSR configurations 
are relevant to the presented study. We can begin with the computation of the path lengths for 
the RSL configuration. The geometry for calculating the Dubins path is shown in Figure 2.11. 


Variable 

Definition 

Units 

(^oAo) 

initial waypoint of search leg 

[none] 


final waypoint of search leg 

[none] 

Yj 

initial heading of searcher 

[radians] 

Yj+I 

final heading of searcher 

[radians] 

r 

minimum turn radius of the searcher 

[meters] 

d 

length of Dubins curve traveled between two waypoints 

[meters] 

D 

configuration of Dubins curve, where D G {RSL^LSR} 

[none] 


distance traveled during initial turn segment 

[meters] 

■'2 

distance traveled during straight-line segment 

[meters] 

■^3 

distance traveled during final turn segment 

[meters] 


pivot point (center) of turn circle for initial turn 

[none] 


pivot point (center) of turn circle for final turn 

[none] 


Table 2.5: Variables associated with the calculation of Dubins curves and their path lengths. 


The searcher is assumed to have the same fixed turn radius for every heading change. This 
results in the straight-line path being along a tangent line between two circles of equal radius. 
There are only four possible tangent lines between any two circles, and because these circles 
are of equal radius, the intersection of the inner tangents occurs at the midpoint between the 
two circle centers; thus, referring to Figure 2.11, we can use the distance between the two pivot 
points, {x^o^,yfo^) and and the turn radius, rtum, to determine the straight-line 

distance and the distance traveled along each curve.The following is a step-by-step derivation of 
the geometry associated with the determination of Dubins distance using a RSL configuration. 
The calculations for the LSR configuration are conducted in an analogous manner. 

We first establish the initial and final turn pivot points, yfo^) and {x^f^,yfj ^), as the offset 
of Tt radians from the heading of the searcher at each point, at a distance of vtum from the Levy 
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Figure 2.11: Illustration showing the variables required in the calculation of the Duhins distance. This 
example shows the geometry for a RSL configuration, hut the process is similar for LSR. 

waypoint (xcTo) or (^/;T/)> shown in Figure 2.12. Having established the pivot points based 
on the Dubins eonfiguration being traveled, we ean now determine the distanee, dist, and the 
offset,0o//, between those pivots points for use in ealeulating the initial turn distance, and 
the straight-line travel distance, 



We next calculate the angle, C, between the intersection point of if and the segment of length 
dist, and the tangent point of if on the first turn circle (^i in Figure 2.11). 



The searcher is assumed to maneuver with a constant turn radius for both the initial and final 
turns, therefore, the angle between tangent point k\ and the pivot point of the final turn is 
always |. Having calculated the value of angles 0o// and C, we can now calculate B, the angle 
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Figure 2.12: Illustration showing the configuration of pivot points at the initial and final waypoints of a 
search leg. The location of the pivot point is at a distance rtum from the waypoint, offset n radians in the 
appropriate direction. 


associated with the arc length of as 

with the distance traveled in the initial turn simply being the length of the arc formed between 
the initial waypoint {xo,yo) and the tangent point ki, or 

sf^ = Brturn. (2.17) 


By a simple application of the Pythagorean theorem we calculate the straight-line distance as 


^2 



(2.18) 


The distance traveled during the second turn depends specifically on the desired final heading 
of the searcher, j8. Again, because rtum is equal for both initial and final turns, symmetry 
exists within the geometry of the Dubins path. In particular the angle, B, and arc length, 
between the tangent point and point perpendicular to the pivot are of equal value in both turns, 
simplifying the calculation of 


sf^ + rturn^, if0<j8<f 

^SL 2;r)], if f < /3 < 


(2.19) 


30 








The distance traveled along the Dubins path is then the sum of the three component parts. 


The calculations for the LSR configuration are conducted in an analogous manner. Note that 
due to the turns occurring in opposite directions from the RSL, the ranges of j8 are different. 
The resulting path length, for the LSR Dubins path is summarized as: 


^LSR 


^LSR 

^2 


.LSR 

^3 



turm 


^ + sm 


-1 1 yrf~yf 

^turn 


+ r[27t-l5], 


if ;r < j8 < ^ 
if ^ < j8 < 2;r 


( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 


The Dubins distance traveled between the waypoints accounts for the constraints on turn radius 
observed in real-world flight, and results in a longer duration of travel than in the instantaneous 
Levy search where no turn radius is implemented. It is necessary to account for this difference 
in distance to ensure accurate calculation of time of flight, ultimately relating to an accurate 
depiction of time to target detection, time to counter-targeting of the searcher, and coverage 
rate. The implementation of the distance calculations in simulation is discussed in Chapter 3. 
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CHAPTER 3: 
Modeling and Simulation 


The flight pattern modeled in this researeh is a eonstant speed trajeetory eomprising segment 
lengths drawn from a Levy distribution. Sueh a model offers an initial mathematieal basis for ex¬ 
ploring the balanee between linear ideal seareh (for improved seareh effeetiveness) and random 
seareh (for improved unpredietability in motion). This ehapter details the funetional proeess of 
the eomputer simulation developed for this researeh. The general flow of the algorithm is repre¬ 
sented by the bloek diagram in Figure 3. land deseribed in pseudo-eode in Algorithm 3.1. The 
assumptions applied during the simulation and a more detailed deseription of eaeh supporting 
funetion follow. 



Figure 3.1: Block diagram showing the relationship and ordering of the various Matlab functions utilized 
throughout execution of the LookAheadLevyFlight.m script. 


3.1 Simulation Model Assumptions 

Although effort is taken to ensure realistie simulation parameters, several assumptions are ap¬ 
plied to the model, distinguishing the simplistie dynamies of this simulation from the employ¬ 
ment of a physieal seareher in a real-world environment. 

The first assumption is that the seareher maintains a eonstant altitude and speed throughout 
exeeution of the flight pattern. The sensor footprint at any diserete time step is based on the 
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altitude of the seareher (e.f. Equation 2.5), and the area eoverage with respeet to time is a 
funetion of seareher speed, as expressed in Equation 2.2. Holding these parameters eonstant at 
values representative of the average speed and altitude exhibited during real-world operations 
simplifies the eoverage ealeulations and faeilitates a dependenee on elapsed flight time only. 

The next assumption is that the sensor utilized by the seareher is mounted on some form of 
gimbal deviee. Most fixed-wing aireraft require banking while turning, whieh may result in 
skewing the shape and plaeement of a fixed sensor footprint. The realistie assumption of a 
gimbal deviee whieh keeps the sensor pointed direetly below the seareher at all times eliminates 
the need for eomplieated ealeulations to determine aetual sensor eoverage shape and plaeement. 

As diseussed in detail in Chapter 2, the seareh grid is divided into one meter by one meter 
diserete eells. These grid eells are small relative to the overall seareh area; thus any grid eell 
that is loeated even partially within the sensor footprint is eonsidered searehed for the purposes 
of this researeh. This simplifieation alleviates the need to perform ealeulations on the exaet 
flight path geometry to determine and traek how mueh of an individual grid eell is searehed 
during eaeh pass of the seareher. 

The seareh area is eireular in geometry and is assumed to eontain no obstruetions whieh pro¬ 
hibit the flight of the seareher at its assigned altitude. In addition to supporting the assumption 
of eonstant altitude, this eliminates the need to model or diseuss eommand and eontrol eonnee- 
tivity issues between the seareher and its eontrolling station, whieh may otherwise arise due to 
hindranees in line of sight. 

Eaeh iteration of the flight pattern begins at the eenter of the eireular seareh area; the transit 
from “home base” to the seareh area and baek are not explieitly aeeounted for in this model. 
The seareh area is assumed to be within reasonable transit distanee from the eontrolling station 
and does not affeet the relative enduranee of the UAV. What is not modeled is the UAV’s suseep- 
tibility to eounter-targeting during transit between the “home base” and the area of operations. 

3.2 Developing the Levy Distribution 

The use of Matlab faeilitates existing aeeess to eommon parametrie distribution funetions, sueh 
as the uniform and exponential distributions. Unfortunately, the Eevy distribution is not among 
the list of eommon distributions provided in Matlab, and thus we ean utilize the eumulative 
distribution funetion of the Eevy, given by Equation 3.1 [45]., to generate an array of seareh 
leg lengths whieh ean be aeeessed during the simulation. Pantaleo et al. diseuss four different 
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methods for parameterizing the Levy distribution for use in simulation, of which we choose the 
transformation method to quantify the cumulative distribution function from the closed form 
Levy defined in Section 1.3 [25,26,45]. 


F(x',jl,c) = f L{l;iJ.,c)dl, 

J —oo 

= erfc 


where erfc{x) is the complementary error function, defined as 


2 /■“ _f2 

erfcix) = / e dt. 

VkJx 


(3.1) 


There are an infinite number of lengths which can be returned from the Levy distribution de¬ 
pending on the accuracy to which we want the lengths represented. For simplicity we choose 
to represent each length to the whole meter (1.0 m). More precise measurements of leg lengths 
were attempted, up to 0.01 meters, but resulted in excessive simulation run times and memory 
use, e.g. the Levy CDF using 1 m precision took 0.28 seconds to build, whereas the 0.01 m pre¬ 
cision CDF took 27.8 seconds, and 0.001 m precision 276.4 seconds. The one meter precision 
coincides with the discretization of the search area and, due to the relative size of the search 
area, provides an accurate representation of coverage ability. We limit the searcher to opera¬ 
tions within the circular search area, and thus a distance greater than twice the radius of that 
area is never traveled on a single straight-line leg, which represents an upper bound on possible 
leg lengths. To address a lower bound on feasible leg lengths, we note that the assumed imple¬ 
mentation of Dubins curves to simulate realistic flight physics requires waypoints separated by 
a minimum of four times the turn radius, although, in general, this minimum distance can be set 
to zero [34]. We set the value of the shift parameter, ji, to this value to represent this minimum 
distance. Although the transform algorithm incorporated here allows for the storage of varying 
degrees of precision on Levy lengths, we calculate only the cumulative probabilities of lengths 
within the defined range bound by the rounded values of /i to 2Rsearch^ to limit the run time and 
memory usage of the simulation when building the Levy CDF. 

We apply each of these values to the Levy CDF by inputting an array, x, of lengths from /i to 
^search, the Scaling parameter, c, and the shift parameter,/!, to the LevyGenerator.m function 
and store the results in a CDF array. Each array index of the Levy CDF maps directly to the 
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index of the original x array of lengths used to build the distribution. To draw a random Levy 
length, a random variable u is drawn from a uniform distribution of zero to one, sueh that 
M = f/(0,1). The Levy CDF array is searehed for this value u, and the index assoeiated with the 
value in the CDF elosest to u is returned, representing the index in the original x array eontaining 
the random Levy length, Ij. This proeess is referred to as inverse transform sampling [46]. 

Figure 3.2 shows the results of this inverse transform sampling for various values of the sealing 
parameter, c, as eompared to the known graph of the Levy PDF. This method provides an 
aeeurate distribution on the Levy lengths for use in the simulation. 



Figure 3.2: Comparison of the known shape of the Levy PDF (a) with the results of the inverse transform 
sampling technique implemented in Matlab for various scaling and shift parameters: (b) /r = 0, c = 0.5, 
(c)n = 0,c = 2.0, (d) At = 0, c = 1.0 


3.3 Restriction of the Searcher to a Finite Area 

The seareher starts eaeh leg of the flight at (Xo^yg), initially assumed the eenter of the seareh area. 
With the length of the seareh leg established, a simple polar distanee ealeulation is eondueted 
to determine the final waypoint {xf,yf) on the eurrent heading a per Equations 2.13 and 2.14. 
If {xf,yf) is within the boundary of the eireular seareh area, the seareher travels along that leg. 
However, if traveling this (Zy, a) eombination results in a violation of the area boundary (i.e., 
outside Rboundary, c.f. Figure 2.7), an adjustment is made to the trajeetory. We propose and 
investigate two methods for eondueting this adjustment. 
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Turn-around Method 

We refer to the first of the two heuristies as the turn-around method. The seareher travels the 
entire length of eaeh randomly drawn Levy leg regardless of the area bounds, whieh may result 
in the seareher temporarily exiting the area. In this adaptive algorithm, the seareher’s position 
with respeet to the area boundary is inspeeted at eaeh waypoint and, upon reeognition of an out- 
of-bounds eondition, direets the seareher baek towards the seareh area via a biased distribution 
on the turn angle, 9. The seareh area is divided into four quadrants equivalent to that of the 
Cartesian plane. Dependent upon whieh quadrant the seareher exeeeds the seareh boundary, the 
limits of the uniform angle distribution are eollapsed to ensure that the next heading plaees the 
seareher on a trajeetory within ±45 degrees of the eenter of the seareh area. This artifieial bias 
on the heading angle is maintained until the seareher is reeognized to be within the bounds of 
the seareh area. A breakdown of the algorithm is shown in Algorithm 3.2. 

By applying the bias to the heading angle, no restrietion is imposed on the lengths of the flight 
legs; thus the distribution of lengths in the flight path remains that of Equation 2.11. Preliminary 
investigation shows that although the ease-dependent heading bias ensures the seareher returns 
to the finite area, the time spent outside of the seareh area eannot be easily predieted. This results 
in large segments of seareh time where no additional area eoverage oeeurs, in effeet wasting the 
seareher’s valuable flight time. Given this operational eonstraint of limited enduranee time in 
praetieal seareh missions, we ehoose to disregard the turn-around method for the remainder of 
this study. 

Look-ahead Method 

We refer to the seeond method for addressing the bounded seareh area as the look-ahead 
method. The seareher again starts at the eenter of the seareh area with random leg lengths 
drawn from Equation 2.11 and random heading ehanges from the uniform distribution on 9. 
Prior to eompleting a seareh leg, the length and angle of the subsequent leg are drawn and a eal- 
eulation eondueted to ensure that transit from the end of the eurrent seareh leg (Ij) to the end of 
the subsequent seareh leg (/;+i) will not plaee the seareher outside of the area boundary. If this 
new length/angle eombination results in a violation of the boundary by the seareher, the length 
and angle are simply disearded and a new length and angle drawn. This proeess is repeated 
until aehieving a length/angle pair whieh meets the eonfinement requirements of the area. This 
proeedure is outlined in Algorithm 3.3. 

The disearding of the lengths that would otherwise plaee the seareher outside the area boundary 
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appears to modify the original probability distribution. The histogram of the leg lengths aetually 
traveled during exeeution of the look-ahead method is shown in Figure 3.3. 



Figure 3.3: The probability histogram of Levy lengths applied to the look-ahead Levy search. The upper 
inlay magnifies the lower values of Ij to show the conformity to the Levy distribution. The lower inlay 
magnifies the lower probability region, showing that the Levy tail is effectively cutoff for values near 
twice the search radius, Rsearch = 6000. 

We plot the CDF of the modified look-ahead Levy distribution versus the standard Levy dis¬ 
tribution of Equation 3.1 and see that indeed, the eontainment proeess results in a more rapid 
rise to maximum probability. The histogram of the look-ahead simulation runs in Figure 3.3(b) 
shows that the general shape of the Levy is maintained, although there is a notieeable inerease 
in probability for the legs less than IRsearch^ thus we enaet the look-ahead eontainment method 
throughout this researeh, with future efforts devoted to examining the effeets of this probability 
shift. The look-ahead method preserves preeious seareh enduranee as compared to the turn¬ 
around method of containment, however the look-ahead is not without fault. The potential 
exists for the searcher to require an infinite number of draws as the length and heading dis¬ 
tributions are discarded due to violation of the boundary conditions. In real world operations 
this could result in inadvertent violation of the boundary by a “runaway searcher” with no up¬ 
dated course vector. As shown in Figure 3.1, we prevent this through implementation of a “stall 
prevention” turn-around function, DubinsTuriiAround.m, after 500 consecutive unsuccessful 
length/angle combinations are drawn. This function simply implements a or 150 degree, 
change in heading angle along the minimum turn radius of the searcher. This value was chosen 
to prevent the searcher becoming stuck in an infinite loop of 180 degree turns, representing a 
“fail-safe” for real-world UAV operations. 
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Figure 3.4: A plot of the CDF compiled from the lengths of the look-ahead Levy (solid) simulation 
versus the standard Levy CDF (dashed) demonstrates the prohahility shift to shorter lengths. 

3.4 Algorithm for Searcher Movement 

3.4.1 Discrete Sensor Sampling 

With the search leg vector determined through the method above, the movement of the searcher 
is implemented in simulation through discrete steps. The searcher position is updated every 
rsensor mctcrs, where Vsensor is the search radius of the sensor. After each position update, the 
LevyDetection.m and LevyCoverage.m functions are run to determine the current operating 
picture. Using this discrete method to determine the amount of area covered accounts for 96% 
of the area that would have been covered through use of a continuous sweep, but saves a signif¬ 
icant amount of simulation run time and processing power, as shown in Figure 3.5. A detailed 
derivation of the calculations used to determine this coverage value can be found in Appendix D. 
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Figure 3.5: This graph shows the effects of varying the spacing between the discretized position updates 
on the simulation runtime. As the spacing increases, the simulation runtime decreases exponentially. 
By choosing a discretized spacing equal to the search radius, we achieve 96% of the coverage attained 
through a continuous search for a fraction of the runtime cost. 
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3.4.2 Movement Along the Flight Path 


The implementation of Dubins curves is an important aspect to modeling real-world flight 
physics throughout this simulation. Once a leg length, Ij, is defined that meets the criteria set 
forth in Section 2.3.2, the initial and final waypoints and headings are passed to the Dubins.m 
function. Calculations are performed within this function to determine the optimal Dubins path 
creating the shortest distance between {xo,yo) and {xf,yf). The mathematical derivation of the 
Dubins calculations is addressed in Section 2.3.3 and the sequence of the Dubins.m function is 
described in Algorithm 3.4. 

DubinsLSR.m and DubinsRSL.m calculate the discretized waypoints on each flight leg. The 
initial and final waypoints, {xo,yo) and {xf,yf), respectively, are passed to the appropriate func¬ 
tion along with the segment lengths of the optimal path configuration as determined in Du¬ 
bins.m. Depending on the path configuration, DubinsLSR.m or DubinsRSL.m generates an 
array of intermediate waypoints along the Dubins path, each with a spacing of rsensor between, 
as discussed in Section 3.4.1. 

Once the waypoint array is returned to the LookAheadLevyFlight.m via Dubins.m and Du¬ 
binsLSR.m or DubinRSL.m, the algorithm begins iterating the array and plotting the way- 
points. The flight time of the searcher is incremented after every intermediate waypoint using 
the distance between those waypoints, AJ and the searcher speed, where 

Ad 

^flight ^flight 3 (3.2) 

As discussed in Section 3.4.1, we have defined the distance to be rsensor- However, we use Ad 
here to account for instances where the distance remaining on a given leg is less than rsensor, as 
accounted for in Algorithm 3.5. 

In Chapter 2, we introduced Tc, or the time the searcher can maneuver in the area before being 
counter-targeted by the adversary. As stated, this time is related to intelligence assessments 
of the technological capability, notably the time the enemy must maintain a target lock on the 
searcher prior to realizing an adequate firing solution. For simplicity, we have assumed that this 
target lock must be achieved during a straight-line path of the searcher; any change in speed 
or heading results in the need to restart targeting efforts. This time to target lock, T^ock, can 
thus be converted into a critical distance, Icrit = Tiqckv, dependent on searcher speed which, if 
exceeded, results in the searcher being counter-targeted. 
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As the algorithm computes the waypoints, the total straight leg distance traveled on the flight 
leg is tracked. If this total straight leg distance exceeds lent, the searcher is considered counter- 
targeted, resulting in a mission failure. The flight time up to that waypoint is recorded as the 
time to counter-targeting, Tc, and the simulation is ended. 

For each of the discretized leg waypoints, the existence of the target within the sensor radius 
is checked and an update to the overall area coverage is performed. Upon reaching the last 
coordinate in the waypoint array, a new length Ij and heading change 6 are drawn and the 
process is repeated until target detection occurs, the searcher is counter-targeted, or the UAV 
reaches maximum flight time. 


3.4.3 Levy Detection 

The LevyDetection.m function is divided into two algorithms, one for conducting a Bayesian 
updated search of the area, the other for conducting the search through employment of a perfect 
cookie cutter sensor. In the case of the Bayesian search, each sensor “look” into a grid cell re¬ 
sults in a Bayesian update as discussed in Section 2.2.2. The sum of the probability mass falling 
within the sensor footprint is updated according to Equations 2.9 and 2.10. The change in prob¬ 
ability resulting from this a posteriori update is dispersed throughout the entire search area, as 
changes in probability proportional to the existing, or a priori, probability in each cell. Based 
on this update, a target detection occurs when the a posteriori probability in a cell exceeds a 
detection confirmation threshold, Over the duration of the search, the concentration 

of probability mass is “pushed” to the vicinity of the target location until, eventually, the prob¬ 
ability of a single grid cell exceeds and the search ends in a successful detection. The 

Bayesian algorithm is also designed to account for sensor imperfection, such as probability of 
missed detection, /3, and probability of false alarm, a. The less perfect the sensor, as indicated 
by higher values of a and j8, the longer the time to target detection tends to be. The summary 
process of the Bayesian update on detection is shown in Algorithm 3.6. 

The second algorithm in the LevyDetection.m function is a perfect detection model to focus 
the analysis on the search pattern itself, shown in Algorithm 3.7. In this algorithm, the target 
location is continually compared to the current searcher location. If the distance between target 
and searcher is less than the sensor range, rsensor, the target is detected. The user specifies the 
type of detection algorithm to be used during the simulation initialization. 
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3.4.4 Levy Loop 

In conjunction with the Bayesian algorithm of the Levy deteetion funetion, the Levy loop fune- 
tion is an optional algorithm intended to deerease the time to target deteetion and eomplete 
area eoverage, whieh is implemented following the LevyDetection.m funetion, as shown in 
Figure 3.1. This funetion is LevyLoop.m divides the updated probability matrix, output from 
the Bayesian algorithm in LevyDetection.m, into equal seetors of size by eells, where N 
is defined by the user during simulation initialization. Eaeh of these seetors represents a lower 
resolution probability mass funetion (eomputationally represented as a matrix) over the entire 
seareh area. 

If no target deteetion oeeurs prior to the expiration of a preset loop timer, whieh begins eounting 
down at the eommeneement of the seareh, a ealeulation is performed to determine whieh of the 
seetors eontains the highest probability mass. This highest probability seetor becomes the target 
destination for the seareher, with heading angles biased to foree the seareher to the area, as 
shown in Algorithm 3.8. When the algorithm determines that the seareher is within the vieinity 
of the high probability area, as determined by the seareher erossing either the horizontal or 
vertieal bounds of the seetor, the loop timer is reset and the seareher resumes the look-ahead 
seareh. 


3.4.5 Levy Coverage 

The LevyCoverage.m funetion takes as an input the eurrent waypoint of the seareher, just as 
in the LevyDetection.m function. The function maintains a Rsearch x R-search coverage matrix, 
where Rsearch is the whole number representation of the search area radius; the discretized 
search area. Initially containing all “zeros”, each matrix block is flagged with a “one” if the 
sensor footprint of the searcher reaches the representative coordinates on the search area. It 
is assumed the sensor footprint is circular; thus only the matrix cells falling within a distance 
r sensor from the searcher position are marked. For each waypoint, the number of ones in the 
matrix is totaled and divided by the total number of cells falling within the A search^ resulting in 
the coverage ratio at the current time of flight, as described in Algorithm 3.4.5. 

, , # matrix cells flagged at tfueht 

F{t) = -- lasAL (33) 

Asearch 
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3.5 Graphical Display 

The simulation is equipped with a graphieal display that may be turned on or off as the user 
desires. The display, a sereen shot of whieh is shown in Figure 3.6, shows the eurrent flight 
trajectory of the searcher, a warning meter to indicate if the searcher is in danger of exceeding 
the critical length, lent, and a plot of coverage ratio versus flight time. The flight trajectory plot 
displays the circular search area as a black line, and the boundary area as a dotted black line. 
The target location (unknown to the searcher) is denoted by an asterisk. The searcher is plotted 
as a dark square surrounded by a lighter colored circle representing the sensor footprint. As the 
searcher moves throughout the area the flight path history and area covered are maintained on 
the plot. 


Levy Search Visualizatior 



Counter-targeting warning 



Distance on leg I. (m) 


Levy Search Coverage Ratio, F(1 



Figure 3.6: Graphical display showing search pattern, coverage ratio, and a warning meter for proximity 
to counter-targeting. 


The Icrit warning meter in the upper right of the display gives a visual representation of the 
probability of the searcher being counter-targeted. As the searcher travels down the straight-leg 
path the meter begins to increase, as represented by a green horizontal bar graph which plots the 
distance traveled on that leg in meters. As the distance traveled exceeds 75% of lent, the meter 
turns yellow. There is no interactive functionality to maneuver the searcher, but it does give the 
user a visual alert that the searcher may be on the verge of counter-targeting. If the straight-leg 
distance exceeds lent, then the meter turns red and the simulation ends. 

The bottom right of the display houses the plot of coverage ratio versus flight time. Each time 
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LevyCoverage.m is run, the computed F{t) is plotted on this graph against the current time of 
flight, in seconds. 

This simulation model encompasses both a visual and mathematical representation of nonde- 
terministic search performance. The graphical output can give the user a feel for the type of 
coverage patterns evolving from such a search, whereas the mathematical output on times to 
target detection and counter-targeting help establish a quantifiable measure on searcher perfor¬ 
mance. In the next chapter we conduct cursory analysis on these measures of performance to 
demonstrate the capabilities inherent to the simulation model. 
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Algorithm 3.1 LookAheadLevyFlight.m 

Build the Levy CDF using LevyGenerator.m 

Initialize seareher position {xo,yo) and heading a, target position, eritieal length, and simu¬ 
lation time 

while Flight time < t^ax do 

Draw Ij from Levy CDF using LevyLength.m and heading ehange 6 from U (0, In) 
Calculate search leg endpoint {xf,yf) and final heading (5 = a +6 
if {xf,yf) does not exceed search area boundary then 

Determine optimal path from {xo,yo) to [xf,yf) using Duhins.m 
while Waypoints remain on current leg do 

if Straight-line leg distance exceeds critical length then 
Searcher counter-targeted. MISSION FAIL 
end if 

if Target located via LevyDetection.m then 
Target detection. MISSION SUCCESS 

end if 

Update F{t) using LevyCoverage.m 

Increment flight time for distance traveled and move to next waypoint on current leg 

end while 
else 

Determine turn around waypoints using DuhinsTurnAround.m 
while Waypoints remain in turn do 

if Target located via LevyDetection.m then 
Target detection. MISSION SUCCESS 

end if 

Update F{t) using LevyCoverage.m 
Move to next waypoint in turn 

Increment flight time for distance traveled and move to next waypoint on current leg 

end while 

a = a + ^ 

end if 

a = p 

Xo =Xf 

yo^yf 

end while 
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Algorithm 3.2 Levy search using the turn-around method 
Initialize searcher position {xo,yo) and heading a 

while t < tmax do 

Draw Ij from Levy CDF using LevyLength.m and heading change 6 
Calculate search leg endpoint {xf,yf) and final heading j8 = a + 0 
Move searcher to (xf^yy) 
if Searcher is outside of area boundary then 

Calculate heading (j) from {xf,yf) to search area center 
Draw Ij from Levy CDF using LevyLength.m 
Draw 6 from biased heading distribution t/(^, |) 

Calculate search leg endpoint {xf^yy) and final heading = ^ + 0 
Move to (;c/,y/) 
end if 

Increment flight time for distance traveled and move to next waypoint on current leg 

end while 


Algorithm 3.3 Levy search using the look-ahead method 
Initialize searcher position {xo,yo) and heading a 

while t < tmax do 

Draw Ij from Levy CDF using LevyLength.m and heading change 0 
Calculate search leg endpoint {xy,yy) and final heading = a + 0 
while4+y}>^Lnd«rydo 

Draw Ij from Levy CDF using LevyLength.m and heading change 0 
Calculate search leg endpoint {xy,yy) and final heading j8 = a + 0 

end while 

Move searcher to {xy,yy) 

Increment flight time for distance traveled and move to next waypoint on current leg 

end while 


Algorithm 3.4 Dubins.m 

Input: lj,{xo,yo),Ci,{xy,yy),li 

ifa<j8<a + 7r then 
Calculate 

Determine the intermediate waypoints using DuhinsRSL.m 
else it a>l5>a — 7t then 
Calculate 

Determine the intermediate waypoints using DuhinsLSR.m 
end if 

Intermediate waypoints returned to LookAheadLevyFlight.m for use in plotting track and 
input to LevyDetection.m and LevyCoverage.m 
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Algorithm 3.5 DubinsLSR.m and DubinsRSL.m, where D G {LSR^RSL} 

Input: if, , if, (;Co, 3^0 ), a, (X/, 3 ;/), j 8 , rsemor 

Initialize segment distance traveled, distance, to zero 
angle = a ± dependent on D 
Determine pivot point of initial turn, (xf,, 3 ^f,) 
while distance < if do 

if if — distance > Vsensor then 

Calculate waypoints along initial turn arc of radius Vtum at a spacing of rsensor 
Store {Xwaypoint ,ywaypoint) in intermediate waypoint array for return to Dubins.m 
angle = angle dependent on D 

distance = distance + rsensor 

else _ 

Calculate final waypoint on initial turn arc at a spacing of if — distance 
Store {xwaypoint ,ywaypoint) in intermediate waypoint array for return to Dubins.m 

, , I s'?—distance , j ^ 7 =^ 

angle = angle — - -, dependent on D 

distance = distance + ^if — distance^ 

end if 
end while 

Initialize segment distance traveled, distance, to zero 
angle = a ± f, dependent on D 
while distance < if do 

if if — distance > Vsensor then 

Calculate waypoints along straight-line segment of heading angle at a spacing of rsensor 
Store {Xsvaypoint,ywaypoint) in intermediate waypoint array 
distance = distance + rsensor 

else 

Calculate final waypoint on straight-line segment at a spacing of if — distance 
Store {Xsvaypoint,ywaypoint) in intermediate waypoint array 

end if 
end while 

Initialize segment distance traveled, distance, to zero 
Determine pivot point of final furn, 

while distance <if do 

if if — distance > Vsensor then 

Calculate waypoinfs along final furn arc of radius rtum at a spacing of rsensor 
Store {x„aypoint,ywaypoint) in intermediate waypoint array 
angle = angle dependent on D 

distance = distance + rsensor 

else 

Calculate final waypoinf on final furn arc af a spacing of i33^ — distance 
Sfore {Xsvaypoint,ywaypoint) in infcrmcdiafe waypoinf array 

angle = j3 

distance = distance + ^if — distance^ 

end if 
end while 




Algorithm 3.6 Bayesian update algorithm incorporated in LevyDetection.m 
Input: prob_density, a, /3, ixpos,ypos),rsensor 
for Each cell c within rsensor of (xpos.ypos) do 

^sensor _„ 

prior ~ l^cPc 

end for 

if (Target < rsensor from {xpos,ypos) AND random t/(0,1) > (no missed detection)) OR 
(random [/(0,1) < a (false detection)) then 

ueiecuon upaaie . pposterior - 

else 


Nondetection update : p^^^^^rior = (T 

end if 

A ^_ -^sensor _ -^sensor 

P posterior P prior 

for All cells c within prob_density do 


posterior (1 - a) (1 )+(j8 )p^ 


j-sensor 
'/I'' prior 
^sensor\\i 
prior 


^ ( _ E. _ 

^posterior ^prior ' ^ prior 


posterior 


prior \ p 


, Vc within r^g, 
, Vc beyond r^, 


of {xpos^yp os) 

.r from {xpos,yp 


Pposterior ^ Pconfirm then 

Target is detected 
To = flight time 

End simulation M/55/OA 5 t/CCE55 

end if 
end for 


Algorithm 3.7 Perfect cookie cutter sensor detection algorithm incorporated in EevyDetec- 

tion.m _ 

Searcher located at position {xpos,ypos) 

Target located at position {xtarget^ytarget) 

if Distance between {xtarget^ytarget) nnd {xpos^ypos) ^sensor then 
Target is detected 
To = flight time 
End simulation 

end if 
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Algorithm 3.8 LevyLoop.m 

Input: sector_prob, a, {xpos,ypos) 

Determine eoordinates of seareh area seetor eontaining highest probability mass, 

{^sector 7 y sector) 

Determine direetion, 7 , of {Xsector,ysector) from eurrent seareher position, {xpos,ypos) 

if 0 < 7 < I then 

while Xpos ^ Xsector OR ypos ^ ysector dO 

Draw Ij from LevyCDF using LevyLength.m 
Draw 6 from biased distribution, t/(0 — 0 + |) 

Move seareher distanee Ij on heading 9 

end while 

else if I < 7 < :;r then 

while Xpos ^ Xsector OR ypos ^ ysector dO 

Draw Ij from LevyCDF using LevyLength.m 
Draw 6 from biased distribution, t/(0 — 0 + |) 

Move seareher distanee Ij on heading 0 

end while 

else if ;r < 7 < ^ then 

while Xpos ^ ^sector OR ypos ^ ysector dO 

Draw Ij from LevyCDF using LevyLength.m 
Draw 0 from biased distribution, t/(0 — 0 + |) 

Move seareher distance Ij on heading 0 

end while 
else 

while Xpos ^ ^sector OR ypos ^ ysector dO 

Draw Ij from LevyCDF using LevyLength.m 
Draw 0 from biased distribution, t/(0 — 0 + |) 

Move searcher distance Ij on heading 0 

end while 
end if 

Reset loop timer to tioop 

Resume search defined in LookAheadLevyFlight.m 
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Algorithm 3.9 LevyCoverage.m 

Input: coverage_plot, {xpos.ypos) ,rsensor 

Coverage eursor {xcoord,ycoord) set at eoordinates (x^aypoint ^sensor-) ywaypoint 
while — Xpos “h sensor do 

ycoord ~ ypos r sensor 
while ycoord ^ ypos “1“ r sensor do 

if y/{Xpos Xcoord)'^ “1“ (ypos ycoord)^) — ^sensor then 

Position {xcoord,ycoord) flagged as searehed on coverage matrix 

end if 

ycoord ycoord “I” ^ 

end while 

^coord ^coord “I” ^ 

end while 

. j # of I’s in coverage matrix 

rea covered — jQjgj number of ceils in search area 
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CHAPTER 4: 

Results, Analysis, and Experimentation 


The purpose of this thesis is the development of a nondeterministie search pattern capable of 
producing area coverage comparable to that of the perfect continuous search, while maintaining 
a sufficient element of randomness conducive to avoiding counter-targeting of the searcher. 
This chapter discusses the results and analysis of simulations performed using various input 
parameters to the Levy distribution and area configuration to develop an understanding of the 
time to detection of a target, captured by the coverage ratio, as well as the time until counter¬ 
targeting, represented as the adversary’s (unknown) capability. 

We first investigate the properties of the search model via an expression for the Levy search 
coverage ratio, which is used to define describe the distribution on the time until the target 
is detected. The inclusion of the counter-targeting element allows for simulation and calcu¬ 
lation of the mission success performance, which is explored through sensitivity analysis on 
the relevant parameters, such as the Levy scaling parameter, c, and the critical time until the 
searcher is counter-targeted, Tiqck- Further discussion of the operational use of these insights 
is developed, presenting capability versus vulnerability curves to be used for decision support. 
To extend this work’s operational relevance, additional study relaxes the assumption of per¬ 
fect sensors by allowing for false positive and false negative detections via a Bayesian search 
framework. Finally, the chapter concludes with a description of live-fly field experiments con¬ 
ducted at Camp Roberts employing the proposed nondeterministie search trajectory generation 
approaches. Implementation using actual tactical UAVs in these actual experiments provide 
proof-of-concept validation of the presented methods. 

4.1 Probability of Mission Success 

Incorporating the primary objective of this research to investigate the impact of the counter¬ 
targeting capability of an adversary against a friendly searching UAV, we define a measure 
of performance referred to as the probability of mission success. Plainly stated, this metric 
quantifies the probability that the searcher finds the stationary target prior to being counter- 
targeted by an enemy weapon system. Recall from Section 2.1 that we refer to the time elapsed 
from the start of the search through target detection as To, and the time elapsed from the start of 
the search through the searcher being counter-targeted as Tc. The probability of mission success 
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is then represented as 


P{Td<Tc). 


(4.1) 


For the purpose of this thesis, we make a distinction between counter-targeting and counterde¬ 
tection. We assume that the territory over which we are searching is under the influence of some 
anti-access/area denial (A2AD) tactic. In other words, the adversary is aware of the searcher’s 
presence in the area (assumed counterdetection), and makes efforts to prevent the searcher from 
locating the target through the use of some anti-air asset (counter-targeting). This asset may 
be a conventional weapon, such as a man-portable air defense system (MANPADS), requiring 
unambiguous travel of the UAV in order for the user to obtain a sufficient target lock and effect 
a hard kill. As countries make technological advances in anti-UAV weaponry, the probability 
exists for increased use of directed energy weapons, requiring a finite amount of laser contact 
time with the surface of the UAV in order to overheat sensitive electronics and effect a mission 
kill. The capability to fire may even be limited by an adversary’s doctrine, requiring a finite 
amount of time or specific number of maneuvers by the targeted UAV before a firing solution 
is considered valid. Recall also from Section 3.4.2, we make the assumption that this finite 
amount of time, annotated Tiqck, is an intelligence assessment of enemy capability, and repre¬ 
sents a random variable with a prescribed probability function. This Tiqck is used to determine 
the maximum straight-line distance, lent, that can be traveled by the searcher prior to being 
counter-targeted. 


4.2 Determining Coverage Ratio 

Recall from Chapter 2 that the coverage ratio is the amount of search area that is uncovered by 
the searcher as a function of time. As discussed in Section 2.1, the continuous perfect search 
results in a linear increase to complete coverage, where as the continuous random search results 
in a logarithmic rise to complete coverage due to the propensity for repeated searches over 
the same grid cell. We compare the performance of the look-ahead Levy search developed in 
this research to these search pattern benchmarks. The discretization of the search area and the 
uniform distribution on the location of the single target allow us to use the coverage ratio as a 
measure of performance for the estimated time to target detection. 

The simulation model is run using the searcher parameters shown in Table 4.1 to determine the 
estimated time required to search the entire area. The results of a single run of the simulation, 
illustrated in Figure 4.1, reveals that in the early stages of the search the coverage rate follows 
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that of the continuous ideal searcher and is nearly linear. Note that there are small increments 
of time which the search algorithm produces results better than what we have described as 
the “perfect” coverage. This is due to imperfections in the search algorithm caused by the 
discretization of the search space, and over time has minimal impact on the search results. As 
the search progresses, the probability of searching grid cells which have already been uncovered 
increases, and the effects of double coverage cause a more logarithmic shape to the curve, 
slowing the coverage rate. 


Variable 

Value 

Units 

V 

25 

[meters/sec] 

c 

0.5 

[none] 


200 

[meters] 

^search 

6000 

[meters] 

^sensor 

161 

[meters] 


Table 4.1: Input parameters to simulation runs used to evaluate coverage rate. 


1 



Figure 4.1: A single run of the simulation model shows that during the early stages of the search, 
a coverage rate similar to that of the continuous ideal searcher is achieved. However, as more area 
is uncovered, the randomness of the pattern causes the searcher to double back over areas previously 
searched, resulting in an logarithmic shape to the curve. 


Given the absence of a closed-form expression for coverage ratio, F{t), the simulation is exe¬ 
cuted n = 1000 times using the same parameters, with the objective of using regression analysis 
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to determine an expression for the (expected) area coverage as a function of time from these 
runs. We conduct a regression analysis on the data to construct a metamodel, denoted F{t), to 
approximate the coverage as a function of time. Noting that the shape of the expected coverate 
ratio mimics a logarithmic rise to complete coverage, similar to that of the random search de¬ 
scribed in Section 2.1, and so we choose the random search equation as the desired parametric 
form. 

To conduct this analysis, we first transform the parametric coverage equation (Equation 2.3) 
into linear form to allow for use of Microsoft Excel’s Data Analysis package. 


Fit) = 1-e"^' 

-log(l-E(t))=At 

We then convert the coverage values obtained from the simulation runs into this form, by taking 
the natural logarithm of one minus the simulated values. A regression is run with this calculated 
log value as the dependent variable, and the simulated time observations as the independent 
variable. The regressor coefficients are shown in the final coverage ratio equation. Equation 4.2. 

^ Y _ g-(0.0000813«+0.035178) ^4 2) 

The statistical results output from Microsoft Excel are summarized in Table 4.2, and the plot 
of the fitted regression values are plotted against the mean simulated coverage values in Eigure 

4.2. 

The coefficients determined from the regression analysis result in a good statistical fit, as indi¬ 
cated by the ;7-value and 7?-squared values in the output file. We note that in the random search 
case, for the given parameter values, the equivalent coefficient (i.e., the coverage rate) is com¬ 
puted to be . = 7.12 X 10 which is of the same order of magnitude as the fitted Eevy 

^search 

coverage rate. However, attempts at further analysis of the regression coefficients does not offer 
insight as to the relation between the coverage ratio and the physical attributes of the searcher or 
search area. Eurther, using this expression for the coverage ratio as an estimate for time to target 
detection does not account for the possibility that the searcher is actively counter-targeted. 
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SUMMARY OUTPUT 


0 

Regression Statistics 
Multiple R 0.999394 

R Square 0.998789 

Adjusted R Square 0.998789 

Standard Error 0.011409 

Observations 4641 


Regression 

df 

1 

SS 

498.0298 

MS 

498.0298 

F 

3825968 

Significance F 

0 




Residual 

4639 

0.603863 

0.00013 






Total 

4640 

498.6336 








Coefficients 

Standard Error 

tStot 

P-value 

Lower 95% 

Upper 95% 

Lower 95.0% 

Upper 95.0% 

Intercept 

0.035178 

0.000336 

104.5482 

0 

0.034519 

0.035838 

0.034519 

0.035838 

Time 

8.13E-05 

4.16E-08 

1956.008 

0 

8.12E-05 

8.14E-05 

8.12E-05 

8.14E-05 


Table 4.2: Regression output from Excel Data Analysis tool. The low p-value and high R-squared values 
indicate the goodness of fit to the coverage rate. 


4.3 Incorporating Counter-Targeting 

In order to facilitate investigation of the impact of the adversary’s ability to target the searcher, 
numerous simulation studies are performed that include variation on parameters relevant to 
both the Levy-based searcher and the enemy’s targeting capability. For the initial studies, we 
consider values of searcher speed, v, and sensor sweep (half) width, rsensor, to be comparable 
to those exhibited in an average mid-sized UAV flying at a 600 meter altitude with a 30° sensor 
field-of-view (c.f. Figure 2.6). Further, the Levy shift parameter, ji, is held constant and dictated 
by the simplifying assumptions of the Dubins algorithm discussed in Section 2.3.3. The target 
is uniformly randomly positioned in the circular search area. 


Variable 

Value 

Units 

Tlock 

90 

[seconds] 

V 

25 

[meters/sec] 

c 

{0.2, 0.5, 1.0, 3.0} 

[none] 


200 

[meters] 

^search 

(3000, 6000, 12000} 

[meters] 

fsensor 

161 

[meters] 


Table 4.3: Input parameters to simulation runs measuring probability of mission success against an 
anti-UAV asset with a Tlqck capability of 90 seconds. 


Each execution of the simulation continues until one of three possible end states is reached: 
either the target is detected, the searcher is counter-targeted, or the endurance of the searcher 
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Figure 4.2: Plot of the mean coverage results from the n = 1000 simulation runs using the parameters in 
Table 4.1. The fitted regression curve is plotted for comparison. 


times out. The first eondition, a target deteetion, is reeorded as a mission sueeess, whereas 
the latter two eonditions, either eounter-targeting of the seareher or exeeeding of the seareher 
enduranee, registers a mission failure. 

The initial numerieal study explores variations in eontrollable parameters, namely those that 
eould potentially be dietated by the UAV eommander, to identify the impaet of the size of the 
seareh area as well as the sealing parameter of the ehosen Levy distribution. In partieular, given 
previous works studying Levy flights on the infinite plane, the effeet of the finite seareh area 
merits investigation. The time neeessary to aehieve target loek is held eonstant at Tiqck = 90 
seeonds, whieh is a representative duration that assumes perfeet intelligenee of the adversary’s 
anti-air targeting eapability. The simulation model is run n = 1,000 times for eaeh of several 
different parameter eombinations, derived from the values shown in Table 4.3. The results of 
these simulations are summarized in Table 4.4. 

Plots of the probability of mission sueeess versus the seareh area size and magnitude of the 
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c 

Rsearch 

[meters] 

E[ro] 

[minutes] 

E[rc] 

[minutes] 

P{Success) 

0.2 

3000 

56.2 

65.5 

0.53 

0.2 

6000 

87.5 

15.1 

0.04 

0.2 

12000 

23.3 

26.4 

0.02 

0.5 

3000 

65.2 

76.3 

0.41 

0.5 

6000 

62.2 

39 

0.03 

0.5 

12000 

N/A 

26.5 

0.00 

1.0 

3000 

34.1 

41.1 

0.40 

1.0 

6000 

3.3 

9.8 

0.05 

1.0 

12000 

16.9 

14.8 

0.01 

3.0 

3000 

26.8 

30.4 

0.29 

3.0 

6000 

4.5 

7.5 

0.07 

3.0 

12000 

9.7 

9.6 

0.01 


Table 4.4: Results of n = 1000 simulation runs against an anti-UAV asset with a T^ock of 90 sec¬ 
onds. Mission success is achieved upon target detection, with mission failure occurring if the searcher 
is counter-targeted. 


scaling parameter are shown in Figure 4.3. These graphs indieate that as the seareh area radius, 
Rsearch, is inereased for a given sealing parameter, c, there is an exponentially detrimental effeet 
on the probability of mission sueeess. In contrast, when Tiqck is known and constant (though, 
in reality, this is a random variable), ehanges of c for a given Rsearch have mueh less effeet on the 
probability of mission sueeess. Based on these results, we foeus the next series of simulations 
on the relationship between Tiqck and Rsearch to identify the sensitivity of mission sueeess 
on these parameters, whieh may provide further insights into the viability of future theoretieal 
eoverage models of Levy walks in finite areas. 

We eonduet n = 100 replieations of eaeh experimental design point, varying values for Rsearch 
and Tiock, with a fixed value of c = 0.5. The remaining seareh parameters are eonsistent with 
those in Table 4.3. The results of these simulations are presented in Table 4.5 and Figure 4.4. 

The results of the simulations foeusing on the relationship between Rsearch and Tiqck again 
indieate that Rsearch has a more dramatie effeet on the probability of mission sueeess. A pairs 
plot of the data is shown in Figure 4.5. The first item to note on this plot is that as the erit- 
ieal length inereases, meaning a longer duration of T^ock, the probability of mission sueeess 
inereases. An inerease in T^ock simply means that the seareher eould fly straighter longer (i.e., 
Icrit is longer), whieh was, for this set of simulations, eontrolled by holding the Levy sealing 
parameter eonstant for a eonsistent probability distribution on the Levy seareh legs. 
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Figure 4.3: (a) Plot showing the effect on mission success of an increase in Rsearch^ given scaling param¬ 
eter, c, and Tlock- (h) Plot showing the effect on mission success of an increase in scaling parameter, c, 
given Rsearch and Tlock- Changes in Rsearch have a greater effect on the prohahility of mission success. 



Search Area Radius, Rsearch {meters) 



(a) 


(b) 


Figure 4.4: (a) Plot of prohahility of mission success versus search area radius for a given Tlock- (h) 
Plot of the prohahility of mission success versus Tlock for a given search radius. 


Further, the probability of mission sueeess deereases in an exponential manner as the radius of 
the seareh area increases. This is partly imposed by the operating conditions of the simulated 
searcher. By confining the searcher to a finite area, the Levy distribution on leg lengths is 
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Tlock 

Rsearch 

[meters] 

Wo] 

[minutes] 

E[rc] 

[minutes] 

P(Success) 

30 

1000 

1.7 

4.6 

0.87 

30 

2000 

3.2 

4.5 

0.35 

30 

3000 

6.0 

3.8 

0.09 

30 

9000 

N/A 

2.8 

0.00 

60 

1000 

2.4 

1.2 

0.98 

60 

2000 

7.5 

16.9 

0.64 

60 

6000 

14.4 

6.5 

0.01 

60 

9000 

1.3 

4.5 

0.01 

90 

1000 

3.4 

0 

1.00 

90 

2000 

14.9 

14.3 

0.90 

90 

3000 

18.4 

17.5 

0.43 

180 

1000 

3.5 

N/A 

1.00 

180 

2000 

18.7 

N/A 

1.00 

180 

3000 

43.2 

56.1 

0.95 

180 

5000 

35.0 

61.1 

0.36 


Table 4.5: Sample of the results of n = 100 simulation runs against an anti-UAV asset with variations 
on Tlock and Rsearch for constant Levy distribution parameters. Mission success is achieved upon target 
detection, with mission failure occurring if the searcher is counter-targeted. Parameters leading to a 
higher probability of mission success than mission failure are bold-faced. 


effectively cut off at twice the area radius, as discussed in Section 2.3.2. As Rsearch is increased 
for a defined distribution on leg lengths (constant c), more of that distribution is “exposed,” 
increasing the overall possibility that a search leg exceeding Icrit is flown. 


4.4 Providing Decision Support 

Given the insights provided by the numerical studies, we can further develop and easily con¬ 
struct reference diagrams that can provide decision support to commanders assessing the search 
capability versus asset vulnerability trade off in conducting associated surveillance missions. 

The simulation developed throughout this research is useful in that it can quantify measures of 
searcher performance given a concise set of operating parameters. We conduct further simula¬ 
tions, separating To and Tc, in order to determine an analytic means for prediction of searcher 
performance through regression of the simulation results. The intent is to separately measure 
searcher performance with respect to target detection time, based on searcher movement and 
sensor capability, and counter-targeting time, based on adversary capability and searcher move¬ 
ment, then realte their statistics to each other to determine an overall probability of mission 
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Figure 4.5: A pairs plot of the data is generated in R to assist in quantifying the relationship between the 
search area size, target lock time, and the prohahility of mission success. 

success. 

4.4.1 Operational Viability of Nondeterministic Search Patterns 

As a preliminary step to first eonvey taetieal relevanee of nondeterministie seareh patterns as 
proposed by this work, we eompare a nominal study of the Levy-based seareh with the expeeted 
time to deteetion performanee for idealized seareh. To generate a working model of the Levy- 
based eoverage, the simulation is exeeuted for n = 10,000 replieations. For these studies, the 
eounter-targeting eapability of the enemy is disabled in the simulation, sueh that the seareh 
progresses until the stationary target is found. This set of simulations employs the same values 
for seareher speed and sensor eapabilities as shown in Table 4.3, with Rsearch the 6000-meter 
level, derived from the maximum area able to be uneovered by the same seareher performing a 
eontinuous ideal (i.e., non-overlapping) seareh. The shift parameter of the Levy distribution, ji. 
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is again set at four times the turn radiusto facilitate the use of the Dubins path algorithms, and 
Levy scaling parameter, c, is set to 0.5. A histogram on the times to target detection is presented 
in Figure 4.6. 



Figure 4.6: Histogram on time to target detection. To, resulting from n = 10,000 simulation runs without 
enemy counter-targeting capability. The dark line function represents an exponential regression fit on 
the distribution. 


We conduct a regression analysis on the distribution of these results, applying an exponential 
fit to model this probability density function on the time to (first) target detection for the given 
search scenario parameters and noting the assumption of the “cookie-cutter” detection model. 
The resulting fitted model for the probability distribution is given by: 

P(Td = t)= e-0-00009095f_ (43) 

This regression fit is represented by the dark line in Figure 4.6. Having constructed a metamodel 
for the probability distribution function for this set of searcher parameters, we can compute the 
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estimated time to target deteetion, E[r£)]. 


E[rz)] 


o 

3 

t e-0.00009095f 


Jo 

3.1 hours. 


(4.4) 


On average, utilizing the parameters established for this set of simulations and ignoring any 
eounter-targeting eapability of the enemy, we expeet the target to be deteeted in 3.1 hours. 
For eomparison, we ealeulate To for a deterministie “spiral-out” seareh pattern illustrated in 
Figure 4.7. This pattern is a variation of the perfeet “lawnmower” seareh pattern deseribed in 
Seetion 2.1, with the seareher beginning at the eenter of the seareh area and eondueting straight- 
line seareh legs marked by 90° turns, moving outwards from the origin until the entire seareh 
area is eleared. The spiraling seareher employed in this eomparison is assigned the same sensor 



Figure 4.7: A variation of the “lawnmower” search discussed in Section 2.1, the “spiral-ouf’ search 
begins at the origin and conducts straight-leg, predictable, non-overlapping sweeps toward the search 
boundary edge. 


and speed eharaeteristies as the seareher in the Levy-based seareh example. In this systematie 
seareh pattern, the time to target deteetion ean be explieitly eomputed (given target loeation) 
using the seareher speed and the sum of the straight-line distanees traveled before the target is 
within range of the sensor. Table 4.6 shows the deterministie time to deteetion for targets plaeed 
at various distanees from the origin. 
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Target distance from origin 
[meters] 

Ideal Tq 
[hours] 

1000 

0.2 

2000 

0.6 

3000 

1.1 

4000 

1.9 

5000 

2.9 

6000 

4.1 


Table 4.6: Times to target detection utilizing a deterministic “spiral-out” search, dependent on distance 
of the target from the beginning of the search. 


The ]E[r£)] of 3.1 hours calculated from the simulation results is within an acceptable magnitude 
for comparison to this deterministic method of search. These results indicate that the detection 
performance of the Levy-based search pattern is comparable to that of the continuous search 
pattern. 


4.4.2 Estimating time to counter-targeting, Tc 

To now provide guidance or “rules of thumb” based on the impact of counter-targeting, we 
can revisit the Levy-based search simulation with varying values of Tiqck- The same searcher 
parameters as previously used (displayed in Table 4.7), are applied to n = 1000 replications 
of the simulation. To specifically investigate Tc, this set of simulations is conducted with no 
target present in the search area, such that the search progresses until a straight-line leg length 
exceeding lent is flown, resulting in the searcher being counter-targeted. 


Variable 

Value 

Units 

Tlock 

{30,60,90,180} 

[seconds] 

V 

25 

[meters/sec] 

c 

0.5 

[none] 


200 

[meters] 

^search 

6000 

[meters] 

fsensor 

161 

[meters] 


Table 4.7: Input parameters for simulation runs to assess Tj) and Tc- 


Figure 4.8 shows the cumulative probability of the searcher being counter-targeted for each 
simulated value of Tiqck, scaled to the upper 20% of the probability density to better display 
the results. This graph represents the probability that counter-targeting of the searcher occurs at 
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or before time t for a given Tiqck- We can use this cumulative distribution graph to determine 
a probability of mission success. 



Search duration, t (hours) 

Figure 4.8: This summary plot of the cumulative distribution functions for varying values of target lock 
time {n = 1000) can be used to determine a measure of the probability of mission success. 

As an example, we plot the expected time to target detection, E[r£)] calculated in Equation 4.4, 
against the cumulative probability of counter-targeting versus an anti-UAV asset with a Tiqck = 
180 seconds. We can read off a surrogate measure of the probability of mission success (c.f. 
Equation 4.1) by using the expected time to detection, which is known a priori by the mission 
commander, since the search parameters are nominally dictated by the search platform: 

P{E[Td] <Tc) = l -P{Tc > E[ro]) 

Using E[r£)] =3.1 hours, we can estimate, based on the expected value graphs of the probability 
of counter-targeting (Eigure 4.8), that the probability of mission success is 1 —0.86 = 0.14 or 
14%. 

Referring back to Table 4.5, for a Tiqck of 180 seconds, the expected probability of mission 
success in a search area of 6000 meter radius is 20%. As expected, the two methods to determine 
this measure of performance produce fairly comparable results, and as such, provide the mission 
commander a means to quantify the probability of mission success for a given set of searcher 
and search area parameters. 
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4.5 Field Experimentation 

In this thesis research, we leveraged unique opportunities to conduct real-world flight testing of 
the Levy-based search algorithm in live-fly contexts. As part of the SECNAV’s initiative for the 
Consortium for Robotics and Unmanned Systems Education and Research (CRUSER) [47] and 
the Naval Postgraduate School’s Advanced Robotic Systems Engineering Eaboratory (AR- 
SENE), live-fly field experimentation was conducted at McMillian Airfield, Camp Roberts, CA. 
The demonstration of the search pattern was conducted as part of a multi-UAV search experi¬ 
ment during the 13-2 Joint Inter-agency Eield Experimentation (HEX) events, an effort by the 
Naval Postgraduate School to explore capability gap solutions and new technologies in support 
of the operational needs homeland defense through discovery, exploration, and experimentation. 

We implement functionality in LookAheadLevyFlight.m to store the major Eevy waypoints 
traversed during the simulated search in Cartesian coordinate form. These “x” and “y” co¬ 
ordinates are input to LevytoLatLon.m, along with a desired starting latitude and longitude 
representing the center of the search area. LevytoLatLon.m iterates the Levy waypoint array 
and converts the distance between each Cartesian waypoint to a Great Circle distance, return¬ 
ing the coinciding latitude and longitude points in an “.fpf” file format for uploading into the 
Procerus Virtual Cockpit (VC) ground control station software. The waypoints in VC are trans¬ 
mitted to the Procerus Unicorn UAV employed during these flight experiments, shown in Eigure 
4.9. 

Eigure 4.10 shows the straight-line path between the major Levy waypoints generated during 
one run of the simulation model. McMillan Airfield is shown to the east of the established 
search area bound within the white dashed circle. The input parameters to each of the four UAV 
flights are listed in Table 4.8. The altitudes of the UAVs were staggered from 300 to 600 meters 
to ensure altitude deconfliction for safety of flight. 


Variable 

Value 

Units 

k^search 

1200 

[meters] 

C 

0.5 

[none] 


200 

[meters] 

V 

15 

[m/s] 


Table 4.8: Searcher parameters established for real-world flight during multi-UAV JIFX demonstration. 

Eigure 4.11 shows the actual flight path taken from the searcher. Appendix C details the pro¬ 
cedure for recovering telemetry data from the Procerus software and converting to a .kml file 
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Figure 4.9: Procerus Unicorn UAV [48] used in multi-UAV look-ahead Levy search demonstration 
conducted at 13-2 JIFX in February, 2013. 

for display in Google Earth. The UAV executes flight trajectories with bounded curvature, con¬ 
sistent with the Dubins curves generated by our simulation, as shown in the magnified flight 
segment displayed in Figure 4.12. However, forcing the UAV to a Dubins flight path requires 
an excessive number of waypoints for any significant duration of flight. The overlay of the 
actual flight trajectory compared to the Dubins path generated by LookAheadLevyFlight.m is 
sufficient evidence that our simulation adequately models the flight physics of a small UAV. 

During the demonstration, ARSENL members were able to employ four UAVs simultaneously 
conducting the look-ahead Eevy search. The flight trajectories are illustrated and distinguished 
by color in Figure 4.13. This experiment represents significant accomplishments in implemen¬ 
tation and system integration, and also highlights the apparent advantage of employing multiple 
searchers evidenced by substantial coverage of the search area. 


4.6 Bayesian Update and Looping Function 

Based on the significant simulation and field experiment studies, key avenues for improving the 
performance of search as well as its operational relevance include the application of Bayesian 
methods to incorporate imperfect sensing and prior intelligence models for likely locations of 
the stationary target. Further, the use of adaptive search, where decisions on where to search are 
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Figure 4.10: Google Earth overlay of major Levy waypoints generated by LookAheadLevyFlight.m 
algorithm. The starting airfield, search center, and search boundary are annotated. 



Figure 4.11: Google Earth overlay of the actual flight path of a single UAV searcher conducting the 
look-ahead Levy search. 


determined based on the results of previous search effort, is known to be significantly advan¬ 
tageous in search [9,49]. We explore such enhancements to the proposed approach, enabling 
preliminary assessment of their value to the nondeterministic search process. 
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Figure 4.12: A magnified section of the flight path traveled hy the UAV compared to the waypoints 
generated through the use of Duhins curves in LookAheadLevy.m. The simulation model created for 
this thesis is demonstrably an adequate representation of the flight physics of a small UAV. 


A Bayesian update seheme is applied to the simulation model, with a “looping” ability whieh 
imparts a bias on the searcher direction based on iteratively updated probabilities of target loca¬ 
tion, thereby potentially decreasing the time to target detection. As described in Section 2.2.2, 
this Bayesian update process utilizes the prior probability that a target exists in the given cell 
and updates that probability based on whether a detection occurs when the searcher “looks” in 
that cell. This Bayesian update process allows for the incorporation and analysis of the charac¬ 
teristics of an imperfect sensor, i.e., the probability of false alarm, a, and probability of missed 
detection, /3. We apply the Bayesian update scheme (where probability mass is shifted in accor¬ 
dance to observations), though we make use of a perfect sensor in order to gain insight into the 
exploitation of the biasing (i.e., looping) that results from the Bayesian update, as discussed in 
detail in Section 2.2.2. For better visual representation and decreased simulation run time, we 
scale down the search parameters to those shown in Table 4.9. 
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Figure 4.13: Google Earth overlay of the actual flight path of a single UAV searcher conducting the 
look-ahead Levy search. 


Variable 

Range 

Units 

Tlock 

90 

[seconds] 

^sensor 

54 

[meters] 

r turn 

25 

[meters] 


100 

[meters] 

C 

0.5 

[none] 

^search 

1500 

[meters] 

a 

0.0 

[none] 

/3 

0.0 

[none] 

hoop 

300 

[seconds] 


Table 4.9: Inputs to the simulation incorporating a Bayesian update and looping function. 


As shown in Figure 4.14, the “binning” funetion of the LevyLoop.m algorithm traeks the prob¬ 
ability density for diserete areas of the seareh grid. This running ealeulation is used to send 
the seareher to clusters of high probability density upon expiration of the user determined loop 
timer. In these simulation runs we employ a five minute loop timer. 

To compare the performance of the Bayesian looped look-ahead search to that of the standard 
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(a) (b) 

Figure 4.14: (a) Heat map showing the distrihution of prohahility of target presence employing a perfect 
sensor Bayesian search over the area, (h) The prohahility density of larger search area sections is used 
to bias the movement of the searcher upon expiration of a loop timer. 


look-ahead Levy seareh used for the majority of this researeh, we randomly draw a target po¬ 
sition (xtarget^ytarget), and run the simulation multiple times for eaeh algorithm. The results of 
average deteetion time and mission sueeess are shown in Table 4.10. 



P(success) 

To.avg 

Favg{To) 

Standard Look-ahead Levy 

1.00 

1221 

0.59 

Looping Look-ahead Levy 

0.90 

480 

0.35 


Table 4.10: Comparison of performance between the standard Look-ahead Levy search and a look-ahead 
Levy search implementing the looping function and Bayesian update. The looping function reduces the 
average time to target detection, but the forced override of the length distribution sometimes leads to 
exceeding the critical length, resulting in counter-targeting of the searcher, and overall mission failure. 


The looping function overrides the random Levy length and heading angle distributions and 
directly and manually steers the searcher toward the area of highest probability density. As 
shown by the slightly lower probability of success, this looping bias may cause the searcher 
to exceed the critical length and thus be counter-targeted. However, the average time to target 
detection when implementing the looping function is almost one-third of the time of the standard 
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look-ahead Levy search. As an added benefit of the looping function, on average, only half as 
much of the search area needed to be uncovered to reveal the target location. This improvement 
may have implications on counter-targeting probability itself, due to less “fly-over” of contested 
territory. These results reiterate the need for accurate intelligence regarding enemy capability. 
By including the looping function into the look-ahead Levy search, but limiting the maximum 
straight-line distance flown by the searcher during the biased relocation to lent, the possibility 
exists to reduce target detection time and exposure to adversarial assets while maximizing the 
probability of mission success. Such assessment and analysis is left for follow-on studies that 
can leverage these benefits. 
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CHAPTER 5: 

Conclusions and Future Avenues 


5.1 Conclusion 

The work accomplished in this thesis provides insight as to the utilization of nondeterministic 
search patterns to counter anti-access area denial (A2AD) tactics becoming more prevalent in 
foreign operations. The Levy distribution on leg lengths serves as only one particular avenue 
to achieving a successful search over this contested territory. The benefits against counter¬ 
targeting of the searcher afforded by the Levy distribution reside in the ability of the user to 
adjust the Levy scaling parameter, c. Recall from Section 1.3 that as the value of c is increased, 
the uniformity of leg lengths returned from the Levy distribution increases, resulting in a longer 
average search leg length. For these higher values of c, the searcher has a greater tendency 
to traverse the entirety of the search area during a span of limited endurance, likely resulting 
in higher percentage of area covered, F{t), for a given search time. However, traveling these 
longer leg lengths results in an increased probability of being counter-targeted due to exceeding 
Icrit- For the sizes of search area explored in this research, changes in c provided only small 
differences in the probability of mission success, albeit with a linear trend indicating that large 
c over a large search area will result in a higher probability of mission failure. A small c 
value (c < 1) provides the highest probability of mission success operating under the parameters 
simulated in this research, since this ensures that the majority of Levy probability density is near 
the /i value, allowing for the implementation of Dubins curves and retaining short flight legs, 
thereby minimizing the probability of exceeding the critical length of the target lock time. 

The Levy distribution has been established as the optimal search pattern when locating mo¬ 
bile targets over an infinitely bound area [18,26]. However, the limitations imposed on the 
searcher’s travel by the finite area boundary used for this thesis may impact the optimality of 
the Levy search, as shown by little difference in performance for changes in Levy parameters. 
Future works focusing on high endurance searchers capable of searching long ranges must be 
undertaken prior to labeling the Levy as the optimal distribution on leg lengths for the scenario 
set here in. Thus, the main contribution of this thesis is the development of the simulation 
model capable of examining various distributions on search leg length and searcher heading. 
The scenario specific to this research focused on an airborne search of a flatland area. However, 
this model has been developed with a robust searcher parameterization capability, and can be 
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applicable to a variety of searcher platforms, from airborne to undersea. 


5.2 Avenues for Future Work 

The development of a simulation capable of analyzing performance for various nondeterministic 
search pattern distributions encompassed the majority of the duration allotted for this research. 

As such, there exists several avenues of future work available for continuation of the efforts 
described in this thesis. An important aspect of a search simulation is the incorporation of 
realism into the model, both operationally and environmentally. There are several extensions 
that can be added to this simulation model to improve upon that realism. 

From an operational perspective, the incorporation of multiple searchers and multiple targets, 
both moving and stationary, will add increased realism and operational relativity. The Levy 
distribution has been shown to outperform the “lawnmower” style continuous search, as well as 
Brownian searches, when attempting to locate multiple moving targets using multiple searchers [23, 
24,34,50,51]. The field experiments conducted for this thesis (Figure 4.13) demonstrate that 
the simultaneous employment of multiple searchers traveling on nondeterministic search paths 
will result in a higher area coverage over a finite time period. This was simply an example of 
multiple non-cooperative searches traveling paths derived from the same algorithm. Introduc¬ 
tion of cooperative functions between these searches could serve to further increase the benefits 
of the nondeterministic search [52]. The incorporation of these cooperative searchers can lead 
to further investigation on proper searcher allocation once a target is found, especially when 
it is expected that multiple targets are within the search area [53,54]. The incorporation and 
analysis of various distributions on searcher leg length and heading angle will provide insight 
as to the optimal nondeterministic pattern when operating in a finite search area. 

Modeling the shape and size of the searcher and targets would greatly effect the product of the 
simulation. Accounting for radar cross section of the searcher relative to a known anti-UAV 
asset can increase fidelity on the measure of counter-targeting time where as increasing the 
size of the target beyond the simple meter squared metric utilized in this research may greatly 
improve time to target detection. 

The enhancement of the operating environment itself will add greater realism to the simulated 
scenario. There have been several studies suggesting algorithmic solutions to obstacle avoid¬ 
ance of unmanned searchers [52,55,56]. Incorporation of physical geography and integration 
of these avoidance algorithms would reveal limitations on UAV flight, and would be particu- 
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larly useful when translating this algorithm into use for a ground searcher. Work accomplished 
regarding optimal coverage using a Dubins vehicle could be extended and incorporated into this 
research by changing the shape of the search area to account for detailed territorial borders or 
geographic constraints establishment, and analyzing the area coverage when a smooth circular 
border is not established [57]. 

Further, adding weather factors such as course drift due to wind or sensor degradation due to 
visibility limitations would also increase the realism of the simulation. 

5.3 Operational Impact and Recommendations 

Since the beginning of the new millennium, the United States has steadily increased reliance 
on unmanned vehicles for conduct of critical missions. These platforms provide a means for 
infiltrating contested territory, for conduct of strike or surveillance operations, without directly 
endangering personnel. However, the world has taken notice of this trend in utilization of un¬ 
manned searchers and, as a result, our adversaries have begin focusing training and technology 
on anti-UAV programs, employing anti-access area denial strategies determined to thwart U.S. 
efforts in intelligence, surveillance, and reconnaissance operations. Perhaps even more signif¬ 
icant, is the subdued sense of conflict associated with destroying these inanimate objects, with 
capturing of U.S. drones being used recently as a bolster of pride, where as the downing and 
capture of a U.S. pilot would most surely be considered an act of war that many countries are 
not willing to instigate. It is imperative to the longevity of these unmanned programs that we 
develop tactics aimed at mitigating the A2AD tactics of our adversaries. The consumers of the 
services provided by these UAVs and other unmanned platforms need to take a step back and 
look toward the future, at the increasing threats and how we can best deploy stealth and random¬ 
ness to mitigate them. The insights contrived from this thesis can help lead the discussion as to 
strategies and technologies we must develop in order to maintain the technical and operational 
edge over our adversaries, not only in the development of further unmanned systems, but in our 
own research and engineering of anti-UAV weapons. 

At the tactical level, this thesis develops a simulation model capable of informing the mission 
commander regarding employment of nondeterministic searchers in a contested territory, an 
increasingly common mission in today’s operational environment. The production of decision 
support tools aimed at shaping the answers to questions such as what specific platform is suited 
for a particular search mission, or how many of that platform need be deployed to provide 
optimal results. Perhaps most importantly is the commanders ability to place a quantifiable 
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estimate on the probability of achieving a mission success given a baseline intelligence report 
regarding enemy capabilities. 
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APPENDIX A: 

Summary of Variables and Symbols 


The MATLAB simulation developed for this thesis ineorporates a robust parameterization ea- 
pability, allowing for the exploration and analysis of multiple seareher types, seareh leg/angle 
distributions, and seareh area eharaeteristies. The simulation inputs whieh ean be set by the user 
are listed in Table A.l. The variables used as outputs to indieate seareh performanee are shown 
in Table A.l. 


Table A.l: Simulation variables set by the user 


Variable 

Description 

Display Options 

plot_ on 

flight, path, on 

1. crit. warning, on 

coverage, plot, on 

plot. 3D 

set to ’ r if any graphieal display is desired 

set to ’ r if flight path traek is to be displayed 

set to ’ r if warning meter for exeeeding eritieal length is to 

be displayed 

set to ’ r if area eoverage plot as a funetion of time is to be 
displayed 

set to ’ r if 3D plot showing relative altitude is desired 

Output Options 

generate. Lat. Lon 

filename. Levy 

filename. Dubin 

set to ’1’ if Cartesian waypoints generated during simula¬ 
tion are to be eonverted to deeimal lat/lon and loaded in 

. f pf file for input to Proeerus Virtual Coekpit 
input filename for general look-ahead Levy waypoints - de¬ 
fault is ‘LevyWaypoints.fpf ’ 

input filename for waypoints foreed to adhere to Dubins 
path - default is ’DubinsWaypoints. fpf’ 

Functional Options 

loop, levy 

set to ’1’ to implement looping funetion, sending seareher 
to highest probability density area after defined t. loop 

time 

Continued on next page 
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Table A.l - continued from previous page 

Variable 

Description 

t_ loop 

set time associated with looping function (how long to 
search before moving to area of highest probability density) 

iteration, count 

set the number of simulation runs to conduct for the current 

set of parameters 

coverage, on 

set to ’1’ if coverage calculations are to be completed (set¬ 
ting to ’0’ turns off coverage if desired to save run time and 
focus purely on detection or counter-targeting parameters) 

target, detect, on 

allows user to turn off the target detection piece in order to 
focus on coverage and/or counter-targeting 

bayesian. on 

set to ’ 1 ’ if Bayesian calculations are to be conducted (set¬ 
ting to ’0’ assumes a perfect sensor and saves time/memory 
on detection calculations) 

counter, target, on 

allows user to turn off the target detection piece in order to 
focus on coverage and/or target detection 

Movement Characteristics of Searcher 

V 

speed of searcher in meters per second 

altitude 

altitude of the searcher corresponds to 2000 ft for UAV in 

meters 

r. turn 

minimum turn radius of the searcher in meters 

UAV. endurance 

maximum endurance of the searcher in meters 

c 

scaling parameter, determines the peak distribution of the 
Levy curve, hence the lengths of each search leg 

Sensor Characteristics of Searcher 

f ov 

field of view of the sensor, in conjunction with altitude, 
determines sensor footprint. Default is 30 degrees as em¬ 
ployed by a typical small UAV 

a 

sensor probability of false detection 

b 

sensor probability of missed detection 

p. confirm 

grid cell probability above which a target is confirmed to be 
located in a given cell, c 


Continued on next page 
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Table A.l - continued from previous page 

Variable 

Description 

p_ deny 

grid cell probability below which a target is confirmed not 
to be located in a given cell, c 

Counter-targeting Parameters 

target, lock 

intelligence-based time required for the anti-UAV asset to 
gain target lock and destroy UAV searcher 

Search Area Parameters 

R_ search 

delta, dist 

radius of the circular search area in meters 

defines the distance between "looks" in meters. In a perfect 
simulation this number would be 0, as the search would be 
continuous. In order to expedite the running of the program 
we use a default step distance of r, returning a minimum 
96% coverage of the area flown over during each search leg 
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Table A.2: Variables which are populated throughout the 
simulation with information useful to the analysis of search 
performance 


Variable 

Description 

Coverage Variables 

coverage_plot 

Matrix representing discretized coverage area. Each cell 
will be ’0’ if not searched, ’ 1’ if searched 

coverage 

Matrix to store percent of area covered at the end of each 
search leg, can store multiple iterations of simulation 

prob_ density 

Matrix representing discretized coverage area. Each cell 
will hold probability of target being located in that cell 

time 

Matrix which stores the time stamp at the end of each search 
leg; Indices map directly to coverage in ‘coverage; matrix 

Performance Metrics 

T_d_array 

Array to store final time to detect the target (‘0’ if no detec¬ 
tion occurs) 

T_c_array 

Array to store final time for searcher to be counter-targeted 
(‘0’ if no counter-targeting occurs) 

Nondeterministic Search Parameters 

angles 

Array to store angles traveled as a result of heading changes 

drawn from the uniform distribution 

distances 

Array to store Eevy lengths applied to each search leg 

DubinsWaypoints 

Matrix to store all v-coordinates, y-coordinates, and altitude 
of intermediate waypoints along Dubins path 

LevyWaypoints 

Matrix to store all v-coordinates, y-coordinates, and altitude 
of major Eevy waypoints only 
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APPENDIX B: 

Software and Hardware 


The following computing software and hardware was used in the development of the search 
model and analysis of simulation results: 

Software: 

• Excel® 

Microsoft®Office®Professional Plus 2010 

• RStudio Statistical Software 
Version 0.95.258 

• Mathworks®MATLAB®R2012a 
Version 7.14.0.739 

The following MATLAB code was generated during development of this simulation. The 
description of each function is discussed in detail in Chapter 3. Each can be made avail¬ 
able upon request at http: //faculty. nps . edu/thchung/ 

■ EookAheadEevyElight.m 

■ EevyGenerator.m 

■ EevyEength.m 

• Dubins.m 

■ DubinsESR.m 

■ DubinsRSE.m 

• DubinsTumAround.m 

■ EevyDetection.m 

■ EevyCoverage.m 

■ EevyEoop.m 

• EevytoEatEon.m 

■ circle.m 

■ spasf.m 
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Hardware: 


• ARSENL-Lab6 

Intel®Core™i7-2700K CPU @ 3.5 GHZ 
32.0 GB 

• Mac-C1MHH11PDV13 
Intel®CoreTMi5 CPU @ 2.4 GHZ 
4.0 GB 

• PC-DELLXPS8300 
Intel®Core™i7-2600K CPU @3.4 GHZ 
12.0 GB 


The following software was used in the conversion of simulation data to real world GPS data for 
the conduct of field experimentation and the display of the resulting telemetry data: Software: 

• Procerus Technologies®Virtual CockpiUM 
Version 2.6.0 build 3 protocol 8.0.1 

• Google®Earth© 

Version 6.2.2.6613 

• GPSVisualizer.com® 

http://www.gpsvisualizer.com 

Hardware: 

• Procerus Unicorn UAV 
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APPENDIX C: 

Waypoint Data Processing for Field Experiments 


This appendix summarizes the proeedures for generating waypoint data for use in the Camp 
Roberts field experiments, as well as the eonversion of telemetry data eolleeted from those 
experiments to .kml files for display in Google Earth. 

C.l Waypoint Generation 

The generation of a waypoint file in .fpf format, for upload direetly to Procerus Virtual Cockpit, 
is a functionality built directly into the LookAheadLevyFlight.m algorithm. In the section of 
the code commented as “output options” there is a Boolean variable named generate_ Lat _ 
Lon. Setting this variable to ‘1’ activates two print scripts, one which writes the main Levy 
waypoints to an .fpf file, the other that writes every waypoint used in the generation of the 
Dubins path flight to a separate .fpf file. It is important to note when attempting to upload these 
waypoints to Virtual Cockpit that the Dubins file can contain in excess of ten times the number 
of waypoints as the Levy waypoint file. 

The names of the output files can be defined by the user by changing the value of the filename_ 
Levy diWA filename_ Dubin string variables in this same section of the LookAheadLevyFlight.m 
code. The default names of the files are “LevyWaypoints.fpf” and “DubinsWaypoints.fpf”. Run 
the LookAheadLevyFlight.m to create and populate theses files, which will be output to the 
main MATLAB directory. 

To upload the waypoints, in Virtual Cockpit click ‘File’, ‘Load Flight Plan’, then select the 
generated .fpf file from the MATLAB directory. Once uploaded, clicking the up arrow at the 
bottom center of the Virtual Cockpit display screen will reveal the list of waypoints, showing 
speed, altitude, latitude, and longitude. Any necessary adjustments can be made here prior to 
transmitting to UAV. 

C.2 Display of Waypoint Data 

The format of the .fpf data file input to Virtual Cockpit requires that each line of waypoint in¬ 
formation be followed by a blank line, effectively double spacing the waypoints. Unfortunately, 
gpsvisualizer. com cannot recognize continuity between the waypoints if the line is skipped, 
and the resultant .kml file will plot each waypoint, but will not connect a track between them. 
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To prevent this, we developed a macro for Excel which simply removes the blank lines between 
waypoints, as well as the extraneous header data required by Virtual Cockpit, but not suitable 
for display in a .kml file. To compress the .csv waypoint file for transformation to a .kml file, 
open the .fpf file generated by LookAheadLevyFlight.m in a text editor. Copy all the waypoint 
information into a .csv file in Excel, then run the macro. 

Once the UAV has conducted the Eevy flight and the telemetry data has been downloaded, the 
data can be sorted and parsed for use in the generation of a .kml file, able to be displayed as 
waypoint tracks in Google Earth. The telemetry information returned from the UAV contains 
108 columns of data, with one row assigned for each telemetry reading (thousands of rows). To 
generate the .kml file, only the basic flight data is require - Heading, Speed, Altitude, Eatitude, 
and Eongitude. The following steps describe the .kml generation process. 

1. Open the telemetry file in MATEAB. The file is named with the UAV identification and 
date, e.g. telemetry_ 2460 _ 2013_ 02_ 14_ 07.m is telemetry data for UAV 2460 flown 
on Eebruary 14, 2013. Run this file to initialize and populate the telemetry.data variable. 

2. Use the ‘dlmwrite’ command in MATEAB to transfer the telemetry data to an Excel .csv 
file. Ensure the numerical precision is set to ‘8’ to so that the full latitude and longitude 
values are transferred. 

dlmwrite{‘filename.csv’,telemetry.data, ‘delimiter’, ‘precision’,8)', 

3. Open the .csv file in Excel and run the first parsing macro by clicking the ‘Developer’ tab, 
‘Macros’, selecting ‘TEEEMETRY_ PARSE’, and clicking ‘Run’. This macro combines 
all time stamp data into and labels each of the columns of telemetry data, placing all 
of the data in a legible, organized format. This will allow the analyst to choose which 
data is relevant to their needs. NOTE: This macro must be run to clean up the telemetry 
data from MATEAB prior to running the .kml generation macro, ‘PARSE_ POR_ KME_ 
GENERATION’. The VBA code for this macro can be made available from https:// 
f acuity.nps.edu/thchung. 

4. To further parse the data down to only those elements required for .kml generation, run the 
.kml parsing function by again clicking the ‘Developer’ tab, ‘Macros’, and selecting the 
‘PARSE_ EOR_ KME_ GENERATION’ macro and clicking ‘Run’. This macro is shown 
in Eisting C.l and can be made available from https ://faculty .nps .edu/thchung. 
This parses the 108 data columns down to the 5 needed to display the flight path in Google 
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Earth. 

5. The easiest way to eonvert the .esv data to a .kml file is to use an external .kml generator 
service, such as provided by http: //GPSvisualizer. com. This website provides free 
.kml conversion, and offers several options on display, ranging from color and shape of 
waypoints to 3D altitude representation. Detailed instructions for use can be found on the 
website. 

6. Once the file is generated, double click the .kml file and the track appears in Google Earth. 


Listing C.l: Parses telemetry data into compressed .esv for input to gpsvisualizer.com 


Sub PARSE_FOR_KML_GENERATION() 

' Parses telemetry data for 

.kml generation for use in Google Earth 

Range ( "A:AG,AM:DB" ) . Select 

Selection . Delete Shift: = 

Range ("Al") . Select 

xlToLeft 

ActiveCell.FormulaRlCl = 

Range ("Bl"). Select 

"Speed" 

ActiveCell.FormulaRlCl = 

Range ("Cl"). Select 

"Altitude " 

ActiveCell.FormulaRlCl = 

Range ("Dl"). Select 

" Heading" 

ActiveCell.FormulaRlCl = 

Range ("El"). Select 

" Latitude " 

ActiveCell.FormulaRlCl = 

Range ("El"). Select 

End Sub 

" Longitude " 
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APPENDIX D; 

Calculations for Discrete Sample Spacing 


To save runtime and proeessing power during simulation runs, we employ a diserete update 
seheme to the seareher movement. Figure D.2 and the subsequent equations provide the math¬ 
ematical background to these discrete updates as introduced in Section 3.4.1. 



Figure D.l: This graph shows the effects of spacing between the discretized looks on the runtime of the 
simulation. As the spacing increases, the simulation runtime decreases exponentially. Choosing a dis¬ 
cretized spacing equal to the search radius provides 96% of the coverage attained through a continuous 
search for a fraction of the runtime cost. 



Figure D.2: This illustration shows the geometry involved in calculating the sensor overlap observed 
during a discretized search with spacing less than twice the sensor radius. 


^lens — '^^wedge '^■^AYZSi (D-1) 

where Aigns represents the overlapping portion of the sensor footprint between two subsequent 
positions of the searcher. A^^edge represents the wedge shaped portion of the overlap outlined in 
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red in Figure D.2 and A/^yzSi represents the triangular portion of the wedge created by points 
YZSi- We subtract Aone time for every two waypoints on the search leg to achieve the 
total new area covered by the searcher on that leg. The remainder of this section focuses on the 
derivation of the equations to conduct these area calculations. 


The area of the wedge is calculated as follows. 
0 = 2cos 


2n, 


r ^sensor 

Awedee ~ ^sensordfsensordO = 

Jo Jo 


fO r3. 


-dO = 


Setting the distance Ad between the “looks” to the radius of the sensor, rsensor, 
equation to 

2 

Awedge ^sensor 

The area of the triangle region YZSi is calculated by 



simplifies the 
(D.2) 


AayZSi — 



1 

4 




with the substitution of d = rsensor yielding 


Aayzs, = (D.3) 

Substitution and simplification of (D.l) reveals the area of double coverage between two subse¬ 
quent waypoints of spacing r. 


Ai 


ens 




or 

Alens = ^-^^diAr^g^sor (D.4) 

Finally, the total area covered (Aeovered) during a single leg of the discretized search can be 
calculated and compared to the total area that would be covered by the perfect continuous 
search on that same leg (Asearch) to determine the percentage of coverage {Pj) achieved through 
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discretization of the searcher movement. 


Pd 


^covered 

^search 

^covered 

^search 


Kr^, 


2/jr 'sensor 


[Kr\ 


2 

sensor 


- 1 . 22844 „,^,) 
1 . 22844 „,„^) 


2ljr 


sensor 


Pd = 0.9566 (D.5) 

Given the benefits of simulation run time and data storage utilizing this discrete movement 
method, and the relative size of the search area compared to the area covered on each search 
leg, a 96% coverage ratio is sufficient for use in this simulation. 
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